/* $Id: Seq_loc.cpp 667223 2023-05-10 18:02:42Z vasilche $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Author:   Author:  Cliff Clausen, Eugene Vasilchenko
 *
 * File Description:
 *   .......
 *
 * Remark:
 *   This code was originally generated by application DATATOOL
 *   using specifications from the ASN data definition file
 *   'seqloc.asn'.
 *
 * ===========================================================================
 */

#include <ncbi_pch.hpp>
#include <serial/enumvalues.hpp>
#include <objects/general/Int_fuzz.hpp>
#include <objects/seqloc/Seq_point.hpp>
#include <objects/seqloc/Seq_interval.hpp>
#include <objects/seqloc/Packed_seqint.hpp>
#include <objects/seqloc/Packed_seqpnt.hpp>
#include <objects/seqloc/Seq_loc_mix.hpp>
#include <objects/seqloc/Seq_loc_equiv.hpp>
#include <objects/seqloc/Seq_bond.hpp>
#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqfeat/Feat_id.hpp>
#include <objects/misc/error_codes.hpp>
#include <util/range_coll.hpp>
#include <objects/seq/seq_id_handle.hpp>
#include <objects/general/Object_id.hpp>
#include <algorithm>


#define NCBI_USE_ERRCODE_X   Objects_SeqLoc

BEGIN_NCBI_SCOPE
BEGIN_objects_SCOPE // namespace ncbi::objects::


// CSeqLocException
const char* CSeqLocException::GetErrCodeString(void) const
{
    switch (GetErrCode()) {
    case eNotSet:       return "eNotSet";
    case eMultipleId:   return "eMultipleId";
    case eUnsupported:  return "eUnsupported";
    case eBadLocation:  return "eBadLocation";
    case eBadIterator:  return "eBadIterator";
    case eIncomatible:  return "eIncomatible";
    case eOutOfRange:   return "eOutOfRange";
    case eOtherError:   return "eOtherError";
    default:            return CException::GetErrCodeString();
    }
}


// constructors
CSeq_loc::CSeq_loc(E_Choice index)
{
    InvalidateCache();
    switch ( index ) {
    case e_Null:
        {
            SetNull();
            break;
        }
    case e_Empty:
        {
            SetEmpty();
            break;
        }
    case e_Whole:
        {
            SetWhole();
            break;
        }
    case e_Int:
        {
            SetInt();
            break;
        }
    case e_Packed_int:
        {
            SetPacked_int();
            break;
        }
    case e_Pnt:
        {
            SetPnt();
            break;
        }
    case e_Packed_pnt:
        {
            SetPacked_pnt();
            break;
        }
    case e_Mix:
        {
            SetMix();
            break;
        }
    case e_Equiv:
        {
            SetEquiv();
            break;
        }
    case e_Bond:
        {
            SetBond();
            break;
        }
    case e_Feat:
        {
            SetFeat();
            break;
        }
    case e_not_set:
    default:
        break;
    }
}


CSeq_loc::CSeq_loc(TId& id, TPoint point, TStrand strand)
{
    InvalidateCache();
    SetPnt(*new CSeq_point(id, point, strand));
}


CSeq_loc::CSeq_loc(TId& id, const TPoints& points, TStrand strand)
{
    InvalidateCache();
    if ( points.size() == 1 ) {
        SetPnt(*new CSeq_point(id, points.front(), strand));
    } else {
        SetPacked_pnt(*new CPacked_seqpnt(id, points, strand));
    }
}


CSeq_loc::CSeq_loc(TId& id, TPoint from, TPoint to, TStrand strand)
{
    InvalidateCache();
    SetInt(*new CSeq_interval(id, from, to, strand));
}


CSeq_loc::CSeq_loc(TId& id, TRanges ranges, TStrand strand)
{
    InvalidateCache();
    if ( ranges.size() == 1 ) {
        SetInt(*new CSeq_interval(id,
            ranges.front().GetFrom(), ranges.front().GetTo(), strand));
    } else {
        SetPacked_int(*new CPacked_seqint(id, ranges, strand));
    }
}


// destructor
CSeq_loc::~CSeq_loc(void)
{
}


inline
void x_Assign(CInt_fuzz& dst, const CInt_fuzz& src)
{
    switch ( src.Which() ) {
    case CInt_fuzz::e_not_set:
        dst.Reset();
        break;
    case CInt_fuzz::e_P_m:
        dst.SetP_m(src.GetP_m());
        break;
    case CInt_fuzz::e_Range:
        dst.SetRange().SetMin(src.GetRange().GetMin());
        dst.SetRange().SetMax(src.GetRange().GetMax());
        break;
    case CInt_fuzz::e_Pct:
        dst.SetPct(src.GetPct());
        break;
    case CInt_fuzz::e_Lim:
        dst.SetLim(src.GetLim());
        break;
    case CInt_fuzz::e_Alt:
        dst.SetAlt() = src.GetAlt();
        break;
    default:
        NCBI_THROW(CSeqLocException, eNotSet,
                   "Int-fuzz is not set");
    }
}


inline
void x_Assign(CSeq_interval& dst, const CSeq_interval& src)
{
    dst.SetFrom(src.GetFrom());
    dst.SetTo(src.GetTo());
    if ( src.IsSetStrand() ) {
        dst.SetStrand(src.GetStrand());
    }
    else {
        dst.ResetStrand();
    }
    dst.SetId().Assign(src.GetId());
    if ( src.IsSetFuzz_from() ) {
        x_Assign(dst.SetFuzz_from(), src.GetFuzz_from());
    }
    else {
        dst.ResetFuzz_from();
    }
    if ( src.IsSetFuzz_to() ) {
        x_Assign(dst.SetFuzz_to(), src.GetFuzz_to());
    }
    else {
        dst.ResetFuzz_to();
    }
}


inline
void x_Assign(CSeq_point& dst, const CSeq_point& src)
{
    dst.SetPoint(src.GetPoint());
    if ( src.IsSetStrand() ) {
        dst.SetStrand(src.GetStrand());
    }
    else {
        dst.ResetStrand();
    }
    dst.SetId().Assign(src.GetId());
    if ( src.IsSetFuzz() ) {
        x_Assign(dst.SetFuzz(), src.GetFuzz());
    }
    else {
        dst.ResetFuzz();
    }
}


inline
void x_Assign(CPacked_seqint& dst, const CPacked_seqint& src)
{
    CPacked_seqint::Tdata& data = dst.Set();
    data.clear();
    ITERATE ( CPacked_seqint::Tdata, i, src.Get() ) {
        data.push_back(CRef<CSeq_interval>(new CSeq_interval));
        x_Assign(*data.back(), **i);
    }
}


inline
void x_Assign(CPacked_seqpnt& dst, const CPacked_seqpnt& src)
{
    if ( src.IsSetStrand() ) {
        dst.SetStrand(src.GetStrand());
    }
    else {
        dst.ResetStrand();
    }
    dst.SetId().Assign(src.GetId());
    if ( src.IsSetFuzz() ) {
        x_Assign(dst.SetFuzz(), src.GetFuzz());
    }
    else {
        dst.ResetFuzz();
    }
    dst.SetPoints() = src.GetPoints();
}


inline
void x_Assign(CSeq_bond& dst, const CSeq_bond& src)
{
    x_Assign(dst.SetA(), src.GetA());
    if ( src.IsSetB() ) {
        x_Assign(dst.SetB(), src.GetB());
    }
    else {
        dst.ResetB();
    }
}


inline
void x_Assign(CSeq_loc_mix& dst, const CSeq_loc_mix& src)
{
    CSeq_loc_mix::Tdata& data = dst.Set();
    data.clear();
    ITERATE ( CSeq_loc_mix::Tdata, i, src.Get() ) {
        data.push_back(CRef<CSeq_loc>(new CSeq_loc));
        data.back()->Assign(**i);
    }
}


inline
void x_Assign(CSeq_loc_equiv& dst, const CSeq_loc_equiv& src)
{
    CSeq_loc_equiv::Tdata& data = dst.Set();
    data.clear();
    ITERATE ( CSeq_loc_equiv::Tdata, i, src.Get() ) {
        data.push_back(CRef<CSeq_loc>(new CSeq_loc));
        data.back()->Assign(**i);
    }
}


void CSeq_loc::Assign(const CSerialObject& obj, ESerialRecursionMode how)
{
    InvalidateCache();
    if ( GetTypeInfo() == obj.GetThisTypeInfo() ) {
        const CSeq_loc& loc = static_cast<const CSeq_loc&>(obj);
        switch ( loc.Which() ) {
        case e_not_set:
            Reset();
            return;
        case e_Null:
            SetNull();
            return;
        case e_Empty:
            SetEmpty().Assign(loc.GetEmpty());
            return;
        case e_Whole:
            SetWhole().Assign(loc.GetWhole());
            return;
        case e_Int:
            x_Assign(SetInt(), loc.GetInt());
            return;
        case e_Pnt:
            x_Assign(SetPnt(), loc.GetPnt());
            return;
        case e_Packed_int:
            x_Assign(SetPacked_int(), loc.GetPacked_int());
            return;
        case e_Packed_pnt:
            x_Assign(SetPacked_pnt(), loc.GetPacked_pnt());
            return;
        case e_Mix:
            x_Assign(SetMix(), loc.GetMix());
            return;
        case e_Equiv:
            x_Assign(SetEquiv(), loc.GetEquiv());
            return;
        case e_Bond:
            x_Assign(SetBond(), loc.GetBond());
            return;
        case e_Feat:
            SetFeat().Assign(loc.GetFeat());
            return;
        }
    }
    CSerialObject::Assign(obj, how);
}


CSeq_loc::TRange CSeq_loc::x_UpdateTotalRange(void) const
{
    TSeqPos range_from = m_TotalRangeCacheFrom.load(memory_order_acquire);
    if ( range_from == TSeqPos(kDirtyCache) ) {
        const CSeq_id* id = 0;
        TRange range = x_CalculateTotalRangeCheckId(id);
        m_IdCache.store(id, memory_order_release);
        m_TotalRangeCacheToOpen.store(range.GetToOpen(), memory_order_relaxed);
        m_TotalRangeCacheFrom.store(range.GetFrom(), memory_order_release);
        return range;
    }
    else {
        TSeqPos range_to_open = m_TotalRangeCacheToOpen.load(memory_order_relaxed);
        return COpenRange<TSeqPos>(range_from, range_to_open);
    }
}


bool CSeq_loc::x_UpdateId(const CSeq_id*& total_id, const CSeq_id* id,
                          bool may_throw) const
{
    if ( total_id == id ) {
        return true;
    }

    if ( !total_id ) {
        total_id = id;
    } else if ( (id  &&  !total_id->Equals(*id)) ) {
        if (may_throw) {
            NCBI_THROW(CSeqLocException, eMultipleId,
                       "CSeq_loc::GetTotalRange() is not defined "
                       "for seq-loc with several different seq-ids");
        } else {
            return false;
        }
    }
    return true;
}



// returns enclosing location range
// the total range is meaningless if there are several seq-ids
// in the location
CSeq_loc::TRange CSeq_loc::x_CalculateTotalRangeCheckId(const CSeq_id*& id) const
{
    TRange total_range;
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
        {
            // Ignore locations without id
            total_range = TRange::GetEmpty();
            break;
        }
    case e_Empty:
        {
            x_UpdateId(id, &GetEmpty());
            total_range = TRange::GetEmpty();
            break;
        }
    case e_Whole:
        {
            x_UpdateId(id, &GetWhole());
            total_range = TRange::GetWhole();
            break;
        }
    case e_Int:
        {
            const CSeq_interval& loc = GetInt();
            x_UpdateId(id, &loc.GetId());
            total_range.Set(loc.GetFrom(), loc.GetTo());
            break;
        }
    case e_Pnt:
        {
            const CSeq_point& pnt = GetPnt();
            x_UpdateId(id, &pnt.GetId());
            TSeqPos pos = pnt.GetPoint();
            total_range.Set(pos, pos);
            break;
        }
    case e_Packed_int:
        {
            // Check ID of each interval
            const CPacked_seqint& ints = GetPacked_int();
            total_range = TRange::GetEmpty();
            ITERATE ( CPacked_seqint::Tdata, ii, ints.Get() ) {
                const CSeq_interval& loc = **ii;
                x_UpdateId(id, &loc.GetId());
                total_range += TRange(loc.GetFrom(), loc.GetTo());
            }
            break;
        }
    case e_Packed_pnt:
        {
            const CPacked_seqpnt& pnts = GetPacked_pnt();
            x_UpdateId(id, &pnts.GetId());
            total_range = TRange::GetEmpty();
            ITERATE( CPacked_seqpnt::TPoints, pi, pnts.GetPoints() ) {
                TSeqPos pos = *pi;
                total_range += TRange(pos, pos);
            }
            break;
        }
    case e_Mix:
        {
            // Check ID of each sub-location.
            const CSeq_loc_mix& mix = GetMix();
            total_range = TRange::GetEmpty();
            ITERATE( CSeq_loc_mix::Tdata, li, mix.Get() ) {
                // instead of Get... to be able to check ID
                total_range += (*li)->x_CalculateTotalRangeCheckId(id);
            }
            break;
        }
    case e_Equiv:
        {
            // Does it make any sense to GetTotalRange() from an equiv?
            const CSeq_loc_equiv::Tdata& equiv = GetEquiv().Get();
            total_range = TRange::GetEmpty();
            ITERATE( CSeq_loc_equiv::Tdata, li, equiv ) {
                total_range += (*li)->x_CalculateTotalRangeCheckId(id);
            }
            break;
        }
    case e_Bond:
        {
            // Does it make any sense to GetTotalRange() from a bond?
            const CSeq_bond& bond = GetBond();
            const CSeq_point& pointA = bond.GetA();
            x_UpdateId(id, &pointA.GetId());
            TSeqPos pos = pointA.GetPoint();
            total_range = TRange(pos, pos);
            if ( bond.IsSetB() ) {
                const CSeq_point& pointB = bond.GetB();
                x_UpdateId(id, &pointB.GetId());
                pos = pointB.GetPoint();
                total_range += TRange(pos, pos);
            }
            break;
        }

    case e_Feat:
    default:
        NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                       "CSeq_loc::CalculateTotalRange(): "
                       "unsupported location type: "<<
                       SelectionName(Which()));
    }

    return total_range;
}


int CSeq_loc::x_CompareSingleId(const CSeq_loc& loc, const CSeq_id* id1,
                                const CSeq_id* id2,
                                TCompareFlags flags) const
{
    if ( !id1 || !id2 ) {
        NCBI_THROW(CSeqLocException, eMultipleId,
                   "CSeq_loc::Compare(): "
                   "cannot compare locations with several different seq-ids");
    }
    if ( int diff = INT_ID_TO(int, id1->CompareOrdered(*id2)) ) {
        // ids are different - order by them
        return diff;
    }

    TSeqPos from1 = GetStart(eExtreme_Positional);
    TSeqPos to1 = GetStop(eExtreme_Positional);
    TSeqPos from2 = loc.GetStart(eExtreme_Positional);
    TSeqPos to2 = loc.GetStop(eExtreme_Positional);

    // (from > to) means circular location.
    // Any circular location is less than (before) non-circular one.
    // If both are circular, compare them regular way.
    bool circular1 = from1 > to1;
    bool circular2 = from2 > to2;
    if ( int diff = circular2 - circular1 ) {
        return diff;
    }

    // smallest left extreme first
    if ( from1 != from2 ) {
        return from1 < from2? -1: 1;
    }
    // longest feature first
    if ( to1 != to2 ) {
        return to1 > to2? -1: 1;
    }
    if (flags & fCompare_Strand) {
        if (!IsSetStrand()) {
            return loc.IsSetStrand() ? -1 : 0;
        }
        else {
            if (!loc.IsSetStrand()) return 1;
        }
        auto s1 = GetStrand();
        auto s2 = loc.GetStrand();
        if (s1 != s2) return s1 < s2 ? -1 : 1;
    }
    return 0;
}

int CSeq_loc::Compare(const CSeq_loc& loc_arg) const
{
    return Compare(loc_arg, fCompare_Default);
}


int CSeq_loc::Compare(const CSeq_loc& loc_arg, TCompareFlags flags) const
{
    // first try fast single-id comparison
    const CSeq_id* id1 = GetId();
    const CSeq_id* id2 = id1 == NULL ? NULL : loc_arg.GetId();
    if (id1 != NULL  &&  id2 != NULL) {
        return x_CompareSingleId(loc_arg, id1, id2, flags);
    }
    // Slow comparison of ranges on each Seq-id separately.
    CSeq_loc_CI iter1(*this,   CSeq_loc_CI::eEmpty_Allow);
    CSeq_loc_CI iter2(loc_arg, CSeq_loc_CI::eEmpty_Allow);
    
    for ( ; iter1 && iter2; ) {
        CRef<CSeq_loc> loc1, loc2;
        for ( int k = 0; k < 2; ++k ) {
            CSeq_loc_CI& iter = k? iter2: iter1;
            CRef<CSeq_loc>& loc = k? loc2: loc1;
            // skip null locations (no Seq-id)
            while ( iter && iter.GetSeq_id().Which() == CSeq_id::e_not_set ) {
                ++iter;
            }
            if ( !iter ) {
                // all the remaining segments were null -> end of location
                loc = null;
                continue;
            }
            const CSeq_id& id = iter.GetSeq_id();
            loc = const_cast<CSeq_loc*>(&*iter.GetRangeAsSeq_loc());
            // skip all sub-locations with the same Seq-id
            while ( ++iter ) {
                if ( !iter.GetSeq_id().Equals(id) ) {
                    if ( iter.GetSeq_id().Which() == CSeq_id::e_not_set ) {
                        // skip null sub-location
                        continue;
                    }
                    else {
                        // different non-null locations (has Seq-id)
                        break;
                    }
                }
                if ( !loc->IsMix() ) {
                    CRef<CSeq_loc> tmp = loc;
                    loc = new CSeq_loc;
                    loc->SetMix().AddSeqLoc(*tmp);
                }
                loc->SetMix().AddSeqLoc(const_cast<CSeq_loc&>(*iter.GetRangeAsSeq_loc()));
            }
        }
        if ( loc1 && !loc2 ) {
            // *this location is longer -> *this > loc
            return 1;
        }
        if ( loc2 && !loc1 ) {
            // second location is longer -> *this < loc
            return -1;
        }
        if ( !loc1 && !loc2 ) {
            // both locations ended
            return 0;
        }
        id1 = loc1->GetId();
        id2 = loc2->GetId();
        if ( int diff = loc1->x_CompareSingleId(*loc2, id1, id2, flags) ) {
            // next sub-locs are different
            return diff;
        }
    }
    if ( iter1 && !iter2 ) {
        // *this location is longer -> *this > loc
        return 1;
    }
    if ( iter2 && !iter1 ) {
        // second locatio is longer -> *this < loc
        return -1;
    }
    // same-length locations
    return 0;
}


void CSeq_loc::PostRead(void) const
{
    InvalidateCache();
}


BEGIN_LOCAL_NAMESPACE;


// Transform a CSeq_loc into a CRange< TSeqPos >,
// that will be its total range.
static inline
CRange<TSeqPos> s_LocToRange(const CSeq_loc& loc )
{
    try {
        return loc.GetTotalRange();
    }
    catch ( CException& /*ignored*/ ) {
        // assume empty
        return CRange<TSeqPos>::GetEmpty();
    }
}


// Transform a CSeq_interval into a CRange< TSeqPos >,
// that will be its total range.
static inline
CRange<TSeqPos> s_LocToRange(const CSeq_interval& interval )
{
    return CRange<TSeqPos>(interval.GetFrom(), interval.GetTo());
};

static inline
bool s_CheckIfMatchesBioseq( const CSeq_loc &loc, 
    const CSeq_loc::ISubLocFilter *filter )
{
    return !filter || (*filter)( loc.GetId() );
}

static inline
bool s_CheckIfMatchesBioseq( const CSeq_interval &loc, 
    const CSeq_loc::ISubLocFilter *filter )
{
    return !filter || (*filter)( &loc.GetId() );
}

// Compare two ranges (reversed on minus strand)
static inline
int s_CompareRanges(const COpenRange<TSeqPos> &x_rng,
                    const COpenRange<TSeqPos> &y_rng,
                    bool minus_strand)
{
    if ( minus_strand ) {
        // This is backwards for compatibliity with C
        // Minus strand features should come in revered order

        // largest right extreme last
        if ( x_rng.GetToOpen() != y_rng.GetToOpen() ) {
            return x_rng.GetToOpen() < y_rng.GetToOpen()? -1: 1;
        }
        // longest last
        if ( x_rng.GetFrom() != y_rng.GetFrom() ) {
            return x_rng.GetFrom() > y_rng.GetFrom()? -1: 1;
        }
    }
    else {
        // smallest left extreme first
        if ( x_rng.GetFrom() != y_rng.GetFrom() ) {
            return x_rng.GetFrom() < y_rng.GetFrom()? -1: 1;
        }
        // longest first
        if ( x_rng.GetToOpen() != y_rng.GetToOpen() ) {
            return x_rng.GetToOpen() > y_rng.GetToOpen()? -1: 1;
        }
    }
    return 0;
}


// Compare two containers with intervals (CSeq_loc or CSeq_interval)
template< class Container1, class Container2 >
int s_CompareIntervals(const Container1& container1,
                       const Container2& container2,
                       bool minus_strand,
                       const CSeq_loc::ISubLocFilter *filter )
{
    typename Container1::const_iterator iter1 = container1.begin();
    typename Container1::const_iterator iter1end = container1.end();
    typename Container2::const_iterator iter2 = container2.begin();
    typename Container2::const_iterator iter2end = container2.end();
    for( ; iter1 != iter1end && iter2 != iter2end ; ++iter1, ++iter2 ) {

        // If specified, skip far location pieces
        if( NULL != filter ) {
            while( iter1 != iter1end && ! s_CheckIfMatchesBioseq(**iter1, filter) ) {
                ++iter1;
            }
            while( iter2 != iter2end && ! s_CheckIfMatchesBioseq(**iter2, filter) ) {
                ++iter2;
            }
            if( iter1 == iter1end || iter2 == iter2end ) {
                break;
            }
        }

        if ( int diff = s_CompareRanges(s_LocToRange(**iter1),
                                        s_LocToRange(**iter2),
                                        minus_strand) ) {
            return diff;
        }
    }

    // finally, shorter sequence first

    if( iter1 == iter1end ) {
        if( iter2 == iter2end ) {
            return 0;
        } else {
            return -1;
        }
    } else {
        return 1;
    }
}


END_LOCAL_NAMESPACE;


int CSeq_loc::CompareSubLoc(const CSeq_loc& loc2, ENa_strand strand, 
    const CSeq_loc::ISubLocFilter *filter) const
{
    bool minus_strand = IsReverse(strand);
    if ( IsMix() ) {
        if ( loc2.IsMix() ) {
            return s_CompareIntervals(GetMix().Get(),
                                      loc2.GetMix().Get(),
                                      minus_strand,
                                      filter );
        }
        else if ( loc2.IsPacked_int() ) {
            return s_CompareIntervals(GetMix().Get(),
                                      loc2.GetPacked_int().Get(),
                                      minus_strand,
                                      filter );
        }
        else {
            // complex loc1 is last on plus strand and first on minus strand
            return minus_strand? -1: 1;
        }
    }
    else if ( IsPacked_int() ) {
        if ( loc2.IsMix() ) {
            return -s_CompareIntervals(loc2.GetMix().Get(),
                                       GetPacked_int().Get(),
                                       minus_strand,
                                       filter );
        }
        else if ( loc2.IsPacked_int() ) {
            return s_CompareIntervals(GetPacked_int().Get(),
                                      loc2.GetPacked_int().Get(),
                                      minus_strand,
                                      filter );
        }
        else {
            // complex loc1 is last on plus strand and first on minus strand
            return minus_strand? -1: 1;
        }
    }
    else {
        if ( loc2.IsMix() || loc2.IsPacked_int() ) {
            // complex loc2 is last on plus strand and first on minus strand
            return minus_strand? 1: -1;
        }
        else {
            // two simple locations
            return 0;
        }
    }
}


bool CSeq_loc::IsSetStrand(EIsSetStrand flag) const
{
    switch ( Which() ) {
    case e_Int:
        return GetInt().IsSetStrand();
    case e_Pnt:
        return GetPnt().IsSetStrand();
    case e_Packed_int:
        return GetPacked_int().IsSetStrand(flag);
    case e_Packed_pnt:
        return GetPacked_pnt().IsSetStrand();
    case e_Mix:
        return GetMix().IsSetStrand(flag);
    case e_Bond:
        return GetBond().IsSetStrand(flag);

    case e_Equiv:
    case e_Feat:
    default:
        return false;
    }
}


ENa_strand CSeq_loc::GetStrand(void) const
{
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
    case e_Empty:
        return eNa_strand_unknown;
    case e_Whole:
        return eNa_strand_both;
    case e_Int:
        return GetInt().IsSetStrand() ? GetInt().GetStrand() : eNa_strand_unknown;
    case e_Pnt:
        return GetPnt().IsSetStrand() ? GetPnt().GetStrand() : eNa_strand_unknown;
    case e_Packed_int:
        return GetPacked_int().GetStrand();
    case e_Packed_pnt:
        return GetPacked_pnt().IsSetStrand() ?
            GetPacked_pnt().GetStrand() : eNa_strand_unknown;
    case e_Mix:
        return GetMix().GetStrand();
    case e_Bond:
        return GetBond().GetStrand();

    case e_Equiv:
    case e_Feat:
    default:
        NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                       "CSeq_loc::GetStrand(): unsupported location type"<<
                       SelectionName(Which()));
    }
}


TSeqPos CSeq_loc::GetStart(ESeqLocExtremes ext) const
{
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
    case e_Empty:
        {
            return kInvalidSeqPos;
        }
    case e_Whole:
        {
            return TRange::GetWhole().GetFrom();
        }
    case e_Int:
        {
            return GetInt().GetStart(ext);
        }
    case e_Pnt:
        {
            return GetPnt().GetPoint();
        }
    case e_Packed_int:
        {
            return GetPacked_int().GetStart(ext);
        }
    case e_Packed_pnt:
        {
            return GetPacked_pnt().GetStart(ext);
        }
    case e_Mix:
        {
            return GetMix().GetStart(ext);
        }
    case e_Bond:
        {
            return GetBond().GetStart(ext);
        }

    case e_Equiv:
    case e_Feat:
    default:
        NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                       "CSeq_loc::GetStart(): "
                       "unsupported location type: "<<SelectionName(Which()));
    }
}


TSeqPos CSeq_loc::GetStop(ESeqLocExtremes ext) const
{
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
    case e_Empty:
        {
            return kInvalidSeqPos;
        }
    case e_Whole:
        {
            return TRange::GetWhole().GetTo();
        }
    case e_Int:
        {
            return GetInt().GetStop(ext);
        }
    case e_Pnt:
        {
            return GetPnt().GetPoint();
        }
    case e_Packed_int:
        {
            return GetPacked_int().GetStop(ext);
        }
    case e_Packed_pnt:
        {
            return GetPacked_pnt().GetStop(ext);
        }
    case e_Mix:
        {
            return GetMix().GetStop(ext);
        }
    case e_Bond:
        {
            return GetBond().GetStop(ext);
        }

    case e_Equiv:
    case e_Feat:
    default:
        NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                       "CSeq_loc::GetStop(): "
                       "unsupported location type: "<<SelectionName(Which()));
    }
}


TSeqPos CSeq_loc::GetCircularLength(TSeqPos seq_len) const
{
    if (seq_len == kInvalidSeqPos) {
        return GetTotalRange().GetLength();
    }
    TSeqPos start = GetStart(eExtreme_Biological);
    TSeqPos stop  = GetStop (eExtreme_Biological);
    bool    minus = IsReverseStrand();

    if (start < stop) {
        return minus ? (seq_len - stop + start + 1) : (stop - start + 1);
    } else {
        return minus ? (start - stop + 1) : (seq_len - start + stop + 1);
    }
}


CSeq_loc::const_iterator CSeq_loc::begin(void) const
{
    return CSeq_loc_CI(*this);
}


CSeq_loc::const_iterator CSeq_loc::end(void) const
{
    return CSeq_loc_CI();
}


// CSeq_loc_CI implementation

SSeq_loc_CI_RangeInfo::SSeq_loc_CI_RangeInfo(void)
    : m_IsSetStrand(false),
      m_Strand(eNa_strand_unknown)
{
}


SSeq_loc_CI_RangeInfo::~SSeq_loc_CI_RangeInfo(void)
{
}


class CSeq_loc_CI_Impl : public CObject
{
public:
    CSeq_loc_CI_Impl(void);
    CSeq_loc_CI_Impl(const CSeq_loc& loc,
                     CSeq_loc_CI::EEmptyFlag empty_flag,
                     CSeq_loc_CI::ESeqLocOrder order);
    virtual ~CSeq_loc_CI_Impl(void) {}

    typedef SSeq_loc_CI_RangeInfo::TRange     TRange;
    typedef vector<SSeq_loc_CI_RangeInfo>     TRanges;

    struct SEquivSet {
        size_t m_StartIndex;

        // parts' end offsets within the equiv set, ordered by position
        // must contain at least one element - the only part
        // all parts must be not empty - contain at least one element
        typedef vector<size_t> TParts;
        typedef TParts::const_iterator TPart_CI;
        typedef TParts::iterator TPart_I;

        TParts m_Parts;

        // return number of Seq-loc element in all equiv parts (always > 0)
        size_t GetElementsCount(void) const
            {
                _ASSERT(!m_Parts.empty());
                _ASSERT(m_Parts.back() > 0);
                return m_Parts.back();
            }
        // return index of the first element in the equiv set
        size_t GetStartIndex(void) const
            {
                return m_StartIndex;
            }
        // return index of element after the equiv set
        size_t GetEndIndex(void) const // exclusive
            {
                return GetStartIndex()+GetElementsCount();
            }

        // return true if the equiv set contains element with specified index
        bool Contains(size_t idx) const
            {
                return idx >= GetStartIndex() && idx < GetEndIndex();
            }

        // return number of equiv parts (always > 0)
        size_t GetPartsCount(void) const
            {
                _ASSERT(!m_Parts.empty());
                return m_Parts.size();
            }

        // return offset of the first element in equiv part within equiv set
        size_t GetPartStartOffset(TPart_CI part) const
            {
                _ASSERT(part < m_Parts.end());
                return part == m_Parts.begin()? 0: *--part;
            }
        // return offset of the end element in equiv part within equiv set
        size_t GetPartEndOffset(TPart_CI part) const
            {
                _ASSERT(part < m_Parts.end());
                return *part;
            }

        // return index of the first element in equiv part within equiv set
        size_t GetPartStartIndex(TPart_CI part) const
            {
                return GetStartIndex()+GetPartStartOffset(part);
            }
        // return index of the end element in equiv part within equiv set
        size_t GetPartEndIndex(TPart_CI part) const
            {
                return GetStartIndex()+GetPartEndOffset(part);
            }

        // return part containing element with specified index
        TPart_I FindPartByElementIndex(size_t idx)
            {
                _ASSERT(Contains(idx));
                return upper_bound(m_Parts.begin(), m_Parts.end(),
                                   idx - GetStartIndex());
            }
        // return part containing element with specified index
        TPart_CI FindPartByElementIndex(size_t idx) const
            {
                _ASSERT(Contains(idx));
                return upper_bound(m_Parts.begin(), m_Parts.end(),
                                   idx - GetStartIndex());
            }
    };
    struct PByLevel {
        bool operator()(const SEquivSet* e1, const SEquivSet* e2) const
            {
                size_t s1 = e1->GetElementsCount();
                size_t s2 = e2->GetElementsCount();
                if ( s1 != s2 ) {
                    // smallest first
                    return s1 < s2;
                }
                s1 = e1->GetPartsCount();
                s2 = e2->GetPartsCount();
                if ( s1 != s2 ) {
                    // enclosing set must have single part with the inner set
                    _ASSERT(s1 == 1 || s2 == 1);
                    // more parts (inner) first
                    return s1 > s2;
                }
                // it's still possible that both equiv sets contain
                // one part with the same set of elements
                // in this case the sets are equivalent and one set contains
                // another, no matter which one, but we still order the sets
                // by their creation order
                return e1 < e2;
            }
    };
    typedef vector<SEquivSet> TEquivSets;

    TRanges& GetRanges(void) { return m_Ranges; }
    const TRanges& GetRanges(void) const { return m_Ranges; }

    bool IsEnd(size_t idx) const { return idx >= m_Ranges.size(); }

    void DeleteRange(size_t idx);
    SSeq_loc_CI_RangeInfo& InsertRange(size_t idx, CSeq_loc::E_Choice type);

    // Point/Interval support
    
    // If any part of the location part is changed the embedding location
    // becomes invalid.
    // This method resets embedding location if it's related to the part only.
    bool WasPoint(const SSeq_loc_CI_RangeInfo& info) const
        {
            return info.m_Loc && info.m_Loc->IsPnt();
        }
    bool ShouldBePoint(const SSeq_loc_CI_RangeInfo& info) const
        {
            if ( !info.m_Loc ) {
                return false;
            }
            switch ( info.m_Loc->Which() ) {
            case CSeq_loc::e_Pnt:
            case CSeq_loc::e_Packed_pnt:
            case CSeq_loc::e_Bond:
                return true;
            default:
                return false;
            }
        }
    // Update CSeq_point object if any
    void UpdatePoint(CSeq_point& pnt,
                     const SSeq_loc_CI_RangeInfo& info) const;
    // Update CSeq_point object if any
    void UpdatePoint(SSeq_loc_CI_RangeInfo& info);
    // Create CSeq_point object
    void SetPoint(SSeq_loc_CI_RangeInfo& info);
    // Update embedded personal Seq-loc object, making sure it's not a point
    void UpdateLoc(SSeq_loc_CI_RangeInfo& info);

    bool CanBePoint(const SSeq_loc_CI_RangeInfo& info) const;

    bool CanBeInterval(const SSeq_loc_CI_RangeInfo& info) const;
    bool CanBePoint(size_t idx) const
        {
            return CanBePoint(m_Ranges[idx]);
        }
    bool CanBeInterval(size_t idx) const
        {
            return CanBeInterval(m_Ranges[idx]);
        }
    bool CanBePacked_pnt(size_t idx_begin, size_t idx_end) const;
    bool CanBePacked_int(size_t idx_begin, size_t idx_end) const;

    CRef<CInt_fuzz> MakeFuzz(const CInt_fuzz& fuzz) const;
    CRef<CSeq_point> MakePoint(const SSeq_loc_CI_RangeInfo& info) const;
    CRef<CSeq_loc> MakeLocPoint(const SSeq_loc_CI_RangeInfo& info) const;

    CRef<CSeq_interval> MakeInterval(const SSeq_loc_CI_RangeInfo& info) const;
    CRef<CSeq_loc> MakeLocInterval(const SSeq_loc_CI_RangeInfo& info) const;
    
    CRef<CSeq_loc> MakeLocPacked_pnt(size_t idx_begin, size_t idx_end) const;
    CRef<CSeq_loc> MakeLocPacked_int(size_t idx_begin, size_t idx_end) const;

    // Bond support
    bool CanBeBond(size_t idx_begin, size_t idx_end) const;
    size_t GetBondBegin(size_t idx) const;
    size_t GetBondEnd(size_t idx) const;
    pair<size_t, size_t> GetBondRange(size_t idx) const
        {
            return make_pair(GetBondBegin(idx), GetBondEnd(idx));
        }
    bool IsInBond(const SSeq_loc_CI_RangeInfo& info) const
        {
            return info.m_Loc && info.m_Loc->IsBond();
        }
    bool IsInBond(size_t idx) const
        {
            return IsInBond(m_Ranges[idx]);
        }
    bool IsBondPartA(size_t idx) const
        {
            return IsInBond(idx) && idx == GetBondBegin(idx);
        }
    bool IsBondPartB(size_t idx) const
        {
            return IsInBond(idx) && idx == GetBondBegin(idx)+1;
        }
    void x_BreakBond(size_t idx)
        {
            SetPoint(m_Ranges[idx]);
        }
    void x_CreateBond(size_t idx)
        {
            CRef<CSeq_loc> loc(new CSeq_loc);
            loc->SetBond();
            m_Ranges[idx].m_Loc = loc;
        }
    void x_AddBondPartB(size_t idx)
        {
            // add part B be to the bond
            m_Ranges[idx+1].m_Loc = m_Ranges[idx].m_Loc;
        }
    void RemoveBond(size_t idx);
    void MakeBondA(size_t idx);
    void MakeBondAB(size_t idx);
    void MakeBondB(size_t idx);
    CRef<CSeq_loc> MakeLocBond(size_t idx_begin, size_t idx_end) const;

    // Equiv support
    bool HasEquivSets(void) const
        {
            return !m_EquivSets.empty();
        }
    bool IsInEquivSet(size_t idx) const
        {
            ITERATE ( TEquivSets, it, m_EquivSets ) {
                if ( it->Contains(idx) ) {
                    return true;
                }
            }
            return false;
        }
    size_t GetEquivSetsCount(size_t idx) const
        {
            size_t count = 0;
            ITERATE ( TEquivSets, it, m_EquivSets ) {
                if ( it->Contains(idx) ) {
                    ++count;
                }
            }
            return count;
        }
    const SEquivSet& GetEquivSet(size_t idx, size_t level) const
        {
            vector<const SEquivSet*> sets;
            ITERATE ( TEquivSets, it, m_EquivSets ) {
                if ( it->Contains(idx) ) {
                    sets.push_back(&*it);
                }
            }
            if ( level >= sets.size() ) {
                NCBI_THROW_FMT(CSeqLocException, eOutOfRange,
                               "CSeq_loc_CI: bad equiv set level: "<<level);
            }
            sort(sets.begin(), sets.end(), PByLevel());
            return *sets[level];
        }
    pair<size_t, size_t> GetEquivSetRange(size_t idx, size_t level) const
        {
            const SEquivSet& s = GetEquivSet(idx, level);
            return make_pair(s.GetStartIndex(), s.GetEndIndex());
        }
    pair<size_t, size_t> GetEquivPartRange(size_t idx, size_t level) const
        {
            const SEquivSet& s = GetEquivSet(idx, level);
            SEquivSet::TPart_CI p = s.FindPartByElementIndex(idx);
            return make_pair(s.GetPartStartIndex(p), s.GetPartEndIndex(p));
        }
    // Return first equiv part bound point that fall into the range.
    // Beginning and ending points are not counted, so 0 is not possible.
    // Return 0 if there are no such breaks.
    size_t HasEquivBreak(size_t begin, size_t end) const;

    typedef set<const SEquivSet*> TUsedEquivs;
    const SEquivSet* FindInnerEquivSet(size_t begin, size_t end,
                                       const TUsedEquivs& used_equivs) const;

    CRef<CSeq_loc> MakeLocOther(const SSeq_loc_CI_RangeInfo& info) const;
    CRef<CSeq_loc> MakeRangeLoc(const SSeq_loc_CI_RangeInfo& info) const;
    CRef<CSeq_loc> MakeLoc(CSeq_loc_I::EMakeType make_type) const;
    CRef<CSeq_loc> MakeLoc(size_t idx_begin,
                           size_t idx_end,
                           CSeq_loc_I::EMakeType make_type,
                           TUsedEquivs& used_equivs) const;

    bool HasChanges(void) const
        {
            return m_HasChanges;
        }

    void SetHasChanges(void)
        {
            m_HasChanges = true;
        }

    CSeq_loc_I::EEquivMode GetEquivMode(void) const
        {
            return m_EquivMode;
        }
    void SetEquivMode(CSeq_loc_I::EEquivMode mode)
        {
            m_EquivMode = mode;
        }

private:
    // Never copy this object
    CSeq_loc_CI_Impl(const CSeq_loc_CI_Impl&);
    CSeq_loc_CI_Impl& operator=(const CSeq_loc_CI_Impl&);

    // Process the location, fill the list
    void x_ProcessLocation(const CSeq_loc& loc);
    void x_ProcessInterval(const CSeq_interval& seq_int,
                           const CSeq_loc& loc);
    void x_ProcessPoint(const CSeq_point& seq_pnt,
                        const CSeq_loc& loc);

    static void x_SetId(SSeq_loc_CI_RangeInfo& info, const CSeq_id& id);

    // Prevent seq-loc destruction
    CConstRef<CSeq_loc>      m_Location;
    // List of intervals
    TRanges                  m_Ranges;
    // List of equiv parts
    TEquivSets               m_EquivSets;
    // Empty locations processing option
    CSeq_loc_CI::EEmptyFlag  m_EmptyFlag;
    // true if anything was changed in the Seq-loc
    bool                     m_HasChanges;
    // equiv editing mode
    CSeq_loc_I::EEquivMode   m_EquivMode;
};


CSeq_loc_CI_Impl::CSeq_loc_CI_Impl(void)
    : m_HasChanges(false),
      m_EquivMode(CSeq_loc_I::eEquiv_none)
{
}


CSeq_loc_CI_Impl::CSeq_loc_CI_Impl(const CSeq_loc&           loc,
                                   CSeq_loc_CI::EEmptyFlag   empty_flag,
                                   CSeq_loc_CI::ESeqLocOrder order)
    : m_Location(&loc),
      m_EmptyFlag(empty_flag),
      m_HasChanges(false),
      m_EquivMode(CSeq_loc_I::eEquiv_none)
{
    x_ProcessLocation(loc);
    if ( order == CSeq_loc_CI::eOrder_Positional  &&  loc.IsReverseStrand() ) {
        reverse(m_Ranges.begin(), m_Ranges.end());
    }
}


void CSeq_loc_CI_Impl::x_SetId(SSeq_loc_CI_RangeInfo& info,
                               const CSeq_id& id)
{
    info.m_Id = &id;
    info.m_IdHandle = CSeq_id_Handle::GetHandle(id);
}


template<class Container>
static void s_ReserveMore(Container& cont, size_t add_size)
{
    if ( cont.size()+add_size > cont.capacity() ) {
        cont.reserve(max(cont.size()+add_size, cont.capacity()*2));
    }
}


void CSeq_loc_CI_Impl::x_ProcessLocation(const CSeq_loc& loc)
{
    switch ( loc.Which() ) {
    case CSeq_loc::e_not_set:
        /*
        NCBI_THROW(CSeqLocException, eNotSet,
                   "CSeq_loc_CI: location is not set");
        */
    case CSeq_loc::e_Null:
    case CSeq_loc::e_Empty:
        {
            if (m_EmptyFlag == CSeq_loc_CI::eEmpty_Allow) {
                SSeq_loc_CI_RangeInfo info;
                if (loc.Which() == CSeq_loc::e_Empty) {
                    x_SetId(info, loc.GetEmpty());
                }
                else {
                    info.m_Id.Reset(new CSeq_id);
                }
                info.m_Range = TRange::GetEmpty();
                info.m_Loc = &loc;
                m_Ranges.push_back(info);
            }
            return;
        }
    case CSeq_loc::e_Whole:
        {
            SSeq_loc_CI_RangeInfo info;
            x_SetId(info, loc.GetWhole());
            info.m_Range = TRange::GetWhole();
            info.m_Loc = &loc;
            m_Ranges.push_back(info);
            return;
        }
    case CSeq_loc::e_Int:
        {
            x_ProcessInterval(loc.GetInt(), loc);
            return;
        }
    case CSeq_loc::e_Pnt:
        {
            x_ProcessPoint(loc.GetPnt(), loc);
            return;
        }
    case CSeq_loc::e_Packed_int:
        {
            const CPacked_seqint::Tdata& data = loc.GetPacked_int().Get();
            s_ReserveMore(m_Ranges, data.size());
            ITERATE ( CPacked_seqint::Tdata, ii, data ) {
                x_ProcessInterval(**ii, loc);
            }
            return;
        }
    case CSeq_loc::e_Packed_pnt:
        {
            const CPacked_seqpnt& pack_pnt = loc.GetPacked_pnt();
            s_ReserveMore(m_Ranges, pack_pnt.GetPoints().size());
            SSeq_loc_CI_RangeInfo info;
            x_SetId(info, pack_pnt.GetId());
            if ( pack_pnt.IsSetStrand() ) {
                info.SetStrand(pack_pnt.GetStrand());
            }
            if (pack_pnt.IsSetFuzz()) {
                info.m_Fuzz.first = info.m_Fuzz.second = &pack_pnt.GetFuzz();
            }
            info.m_Loc = &loc;
            ITERATE ( CPacked_seqpnt::TPoints, pi, pack_pnt.GetPoints() ) {
                info.m_Range.Set(*pi, *pi);
                m_Ranges.push_back(info);
            }
            return;
        }
    case CSeq_loc::e_Mix:
        {
            const CSeq_loc_mix::Tdata& data = loc.GetMix().Get();
            s_ReserveMore(m_Ranges, data.size());
            ITERATE(CSeq_loc_mix::Tdata, li, data) {
                x_ProcessLocation(**li);
            }
            return;
        }
    case CSeq_loc::e_Equiv:
        {
            const CSeq_loc_equiv::Tdata& data = loc.GetEquiv().Get();
            s_ReserveMore(m_Ranges, data.size());
            SEquivSet eq_set;
            size_t start = m_Ranges.size();
            eq_set.m_StartIndex = start;
            ITERATE(CSeq_loc_equiv::Tdata, li, data) {
                size_t part_start = m_Ranges.size();
                x_ProcessLocation(**li);
                size_t part_end = m_Ranges.size();
                if ( part_start != part_end ) {
                    // add only non-empty parts
                    eq_set.m_Parts.push_back(part_end-start);
                }
            }
            if ( !eq_set.m_Parts.empty() ) {
                // add only non-empty equiv sets
                m_EquivSets.push_back(eq_set);
            }
            return;
        }
    case CSeq_loc::e_Bond:
        {
            const CSeq_bond& bond = loc.GetBond();
            x_ProcessPoint(bond.GetA(), loc);
            if ( bond.IsSetB() ) {
                x_ProcessPoint(bond.GetB(), loc);
            }
            return;
        }
    default:
        NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                       "CSeq_loc_CI: unsupported location type: "<<
                       loc.SelectionName(loc.Which()));
    }
}


void CSeq_loc_CI_Impl::x_ProcessInterval(const CSeq_interval& seq_int,
                                         const CSeq_loc& loc)
{
    _ASSERT(loc.Which() == CSeq_loc::e_Int ||
            loc.Which() == CSeq_loc::e_Packed_int);
    SSeq_loc_CI_RangeInfo info;
    x_SetId(info, seq_int.GetId());
    info.m_Range.Set(seq_int.GetFrom(), seq_int.GetTo());
    if ( seq_int.IsSetStrand() ) {
        info.SetStrand(seq_int.GetStrand());
    }
    info.m_Loc = &loc;
    if (seq_int.IsSetFuzz_from()) {
        info.m_Fuzz.first = &seq_int.GetFuzz_from();
    }
    if (seq_int.IsSetFuzz_to()) {
        info.m_Fuzz.second = &seq_int.GetFuzz_to();
    }
    m_Ranges.push_back(info);
}


void CSeq_loc_CI_Impl::x_ProcessPoint(const CSeq_point& seq_pnt,
                                      const CSeq_loc& loc)
{
    _ASSERT(loc.Which() == CSeq_loc::e_Pnt ||
            loc.Which() == CSeq_loc::e_Bond);
    SSeq_loc_CI_RangeInfo info;
    x_SetId(info, seq_pnt.GetId());
    info.m_Range.Set(seq_pnt.GetPoint(), seq_pnt.GetPoint());
    if ( seq_pnt.IsSetStrand() ) {
        info.SetStrand(seq_pnt.GetStrand());
    }
    info.m_Loc = &loc;
    if (seq_pnt.IsSetFuzz()) {
        info.m_Fuzz.first = info.m_Fuzz.second = &seq_pnt.GetFuzz();
    }
    m_Ranges.push_back(info);
}


size_t CSeq_loc_CI_Impl::HasEquivBreak(size_t begin,
                                       size_t end) const
{
    size_t ret = end;
    ITERATE ( TEquivSets, it, m_EquivSets ) {
        const SEquivSet& eq_set = *it;
        size_t eq_begin = eq_set.GetStartIndex();
        size_t eq_end = eq_set.GetEndIndex();
        if ( eq_end <= begin || eq_begin >= end ) {
            // no overlaps
            continue;
        }
        if ( eq_begin > begin && eq_begin < end ) {
            // range overlaps with the beginning of the equiv set
            ret = min(ret, eq_begin);
            continue;
        }
        // otherwise the equiv starts before the range and ends somewhere
        // after beginning of the range.
        // so the first part that contains range's begin point is the one,
        // end its end point is one of the breaking points we are looking for.
        SEquivSet::TPart_CI part = eq_set.FindPartByElementIndex(begin);
        ret = min(ret, eq_set.GetPartEndIndex(part));
    }
    // the first break point may be after the requested range
    return ret == end? 0: ret;
}


const CSeq_loc_CI_Impl::SEquivSet*
CSeq_loc_CI_Impl::FindInnerEquivSet(size_t begin,
                                    size_t end,
                                    const TUsedEquivs& used_equivs) const
{
    const SEquivSet* ret = 0;
    ITERATE ( TEquivSets, it, m_EquivSets ) {
        const SEquivSet& eq_set = *it;
        size_t eq_begin = eq_set.GetStartIndex();
        size_t eq_end = eq_set.GetEndIndex();
        if ( eq_begin < begin || eq_end > end ) {
            continue;
        }
        if ( used_equivs.find(&eq_set) != used_equivs.end() ) {
            // already processed
            continue;
        }
        if ( !ret || PByLevel()(ret, &eq_set) ) {
            ret = &eq_set;
        }
    }
    return ret;
}


void CSeq_loc_CI_Impl::DeleteRange(size_t idx)
{
    SetHasChanges();
    m_Ranges.erase(m_Ranges.begin()+idx);
    // update equiv sets
    ERASE_ITERATE ( TEquivSets, it, m_EquivSets ) {
        SEquivSet& eq_set = *it;
        size_t start_index = eq_set.GetStartIndex();
        if ( idx < start_index ) {
            // shift whole equiv set because an element before is deleted
            --eq_set.m_StartIndex;
            continue;
        }
        size_t offset = idx - start_index;
        size_t prev_offset = 0;
        ERASE_ITERATE ( SEquivSet::TParts, it2, eq_set.m_Parts ) {
            if ( offset < *it2 ) {
                --*it2;
                if ( *it2 == prev_offset ) {
                    VECTOR_ERASE(it2, eq_set.m_Parts);
                }
                else {
                    prev_offset = *it2;
                }
            }
        }
        if ( eq_set.GetElementsCount() == 0 ) {
            VECTOR_ERASE(it, m_EquivSets);
        }
    }
}


static inline bool s_InsertExpands(size_t idx, size_t start, size_t end)
{
    return (idx > start && idx <= end) || (idx == start && idx == 0);
}


SSeq_loc_CI_RangeInfo& CSeq_loc_CI_Impl::InsertRange(size_t idx,
                                                     CSeq_loc::E_Choice type)
{
    SetHasChanges();
    SSeq_loc_CI_RangeInfo& info =
        *m_Ranges.insert(m_Ranges.begin()+idx, SSeq_loc_CI_RangeInfo());
    CConstRef<CSeq_loc> loc;
    // connect with previous or next part (if any)
    if ( idx > 0 ) {
        loc = m_Ranges[idx-1].m_Loc;
    }
    else if ( idx < m_Ranges.size() ) {
        loc = m_Ranges[idx].m_Loc;
    }
    if ( loc ) {
        // if type is incompatible with part it's connected to, do not connect
        switch ( loc->Which() ) {
        case CSeq_loc::e_Null:
        case CSeq_loc::e_Empty:
        case CSeq_loc::e_Whole:
        case CSeq_loc::e_Pnt:
        case CSeq_loc::e_Int:
            loc = null;
            break;
        case CSeq_loc::e_Bond:
            if ( type != CSeq_loc::e_Pnt ) {
                loc = null;
            }
            break;
        case CSeq_loc::e_Packed_pnt:
            if ( type != CSeq_loc::e_Pnt ) {
                loc = null;
            }
            break;
        case CSeq_loc::e_Packed_int:
            if ( type != CSeq_loc::e_Pnt && type != CSeq_loc::e_Int ) {
                loc = null;
            }
            break;
        default:
            break;
        }
    }
    info.m_Loc = loc;

    SEquivSet* add_to_set_beg = 0;
    vector<SEquivSet*> add_to_set_mid;
    SEquivSet* add_to_set_end = 0;
    // update existing equiv sets
    NON_CONST_ITERATE ( TEquivSets, it, m_EquivSets ) {
        SEquivSet& eq_set = *it;
        size_t eq_start = eq_set.GetStartIndex();
        size_t eq_end = eq_set.GetEndIndex();
        if ( idx == eq_start &&
             (GetEquivMode() == CSeq_loc_I::eEquiv_prepend ||
              GetEquivMode() == CSeq_loc_I::eEquiv_new_part) ) {
            // new element is inserted just before the equiv and will be added
            // shift it for now
            ++eq_set.m_StartIndex;
            if ( !add_to_set_beg || PByLevel()(&eq_set, add_to_set_beg) ) {
                add_to_set_beg = &eq_set;
            }
        }
        else if ( idx == eq_end &&
                  (GetEquivMode() == CSeq_loc_I::eEquiv_append ||
                   GetEquivMode() == CSeq_loc_I::eEquiv_new_part) ) {
            // new element is inserted just after the equiv and will be added
            if ( !add_to_set_end || PByLevel()(&eq_set, add_to_set_end) ) {
                add_to_set_end = &eq_set;
            }
        }
        else if ( idx > eq_start && idx < eq_end ) {
            // new element is completely within eq_set, so it will be added
            add_to_set_mid.push_back(&eq_set);
        }
        else if ( idx <= eq_start ) {
            // new element is completely before the equiv - shift right
            ++eq_set.m_StartIndex;
        }
    }
    sort(add_to_set_mid.begin(), add_to_set_mid.end(), PByLevel());
    if ( add_to_set_beg ) {
        // add at the beginning of the set
        NON_CONST_ITERATE ( SEquivSet::TParts, it, add_to_set_beg->m_Parts ) {
            ++*it;
        }
        if ( GetEquivMode() == CSeq_loc_I::eEquiv_new_part ) {
            // add as a new part
            add_to_set_beg->m_Parts.insert(add_to_set_beg->m_Parts.begin(), 1);
            SetEquivMode(CSeq_loc_I::eEquiv_append);
        }
    }
    else if ( add_to_set_end ) {
        // add at the end of the set
        if ( GetEquivMode() == CSeq_loc_I::eEquiv_new_part ) {
            // add as a new part
            add_to_set_end->m_Parts.push_back(add_to_set_end->m_Parts.back()+1);
            SetEquivMode(CSeq_loc_I::eEquiv_append);
        }
        else {
            ++add_to_set_end->m_Parts.back();
        }
    }
    if ( GetEquivMode() == CSeq_loc_I::eEquiv_new_part &&
         !add_to_set_mid.empty() ) {
        SEquivSet& eq_set = **add_to_set_mid.begin();
        add_to_set_mid.erase(add_to_set_mid.begin());
        SEquivSet::TPart_I it = eq_set.FindPartByElementIndex(idx);
        it = eq_set.m_Parts.insert(it, idx-eq_set.GetStartIndex()+1)+1;
        for ( ; it != eq_set.m_Parts.end(); ++it ) {
            ++*it;
        }
        SetEquivMode(CSeq_loc_I::eEquiv_append);
    }
    ITERATE ( vector<SEquivSet*>, set_it, add_to_set_mid ) {
        SEquivSet& eq_set = **set_it;
        SEquivSet::TPart_I it = eq_set.FindPartByElementIndex(idx);
        for ( ; it != eq_set.m_Parts.end(); ++it ) {
            ++*it;
        }
    }
    // now all existing equivs were updated

    if ( GetEquivMode() == CSeq_loc_I::eEquiv_new_equiv ) {
        // create new equiv
        SEquivSet eq_set;
        eq_set.m_StartIndex = idx;
        eq_set.m_Parts.push_back(1);
        m_EquivSets.push_back(eq_set);
        SetEquivMode(CSeq_loc_I::eEquiv_append);
    }
    return info;
}


bool CSeq_loc_CI_Impl::CanBePoint(const SSeq_loc_CI_RangeInfo& info) const
{
    if ( info.m_Range.GetLength() != 1 ) {
        return false;
    }
    if ( info.m_Fuzz.first != info.m_Fuzz.second ) {
        return false;
    }
    if ( !info.m_IdHandle ) {
        return false;
    }
    return true;
}


bool CSeq_loc_CI_Impl::CanBeInterval(const SSeq_loc_CI_RangeInfo& info) const
{
    if ( info.m_Range.Empty() || info.m_Range.IsWhole() ) {
        return false;
    }
    if ( !info.m_IdHandle ) {
        return false;
    }
    return true;
}


bool CSeq_loc_CI_Impl::CanBeBond(size_t idx_begin, size_t idx_end) const
{
    size_t count = idx_end - idx_begin;
    if ( count != 1 && count != 2 ) {
        return false;
    }
    if ( !IsInBond(idx_begin) ) {
        return false;
    }
    if ( GetBondRange(idx_begin) != make_pair(idx_begin, idx_end) ) {
        return false;
    }
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        _ASSERT(IsInBond(idx));
        const SSeq_loc_CI_RangeInfo& info = m_Ranges[idx];
        if ( !CanBePoint(info) ) {
            return false;
        }
    }
    if ( HasEquivBreak(idx_begin, idx_end) ) {
        return false;
    }
    return true;
}


size_t CSeq_loc_CI_Impl::GetBondBegin(size_t idx) const
{
    const CSeq_loc* loc = m_Ranges[idx].m_Loc.GetPointerOrNull();
    _ASSERT(loc && loc->IsBond());
    while ( idx > 0 && m_Ranges[idx-1].m_Loc == loc ) {
        --idx;
    }
    return idx;
}


size_t CSeq_loc_CI_Impl::GetBondEnd(size_t idx) const
{
    const CSeq_loc* loc = m_Ranges[idx].m_Loc.GetPointerOrNull();
    _ASSERT(loc && loc->IsBond());
    while ( idx < m_Ranges.size() && m_Ranges[idx].m_Loc == loc ) {
        ++idx;
    }
    return idx;
}


void CSeq_loc_CI_Impl::RemoveBond(size_t idx)
{
    if ( !IsInBond(idx) ) {
        NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                       "CSeq_loc_I::RemoveBond(): "
                       "there is no bond at current position");
    }
    pair<size_t, size_t> range = GetBondRange(idx);
    SetHasChanges();
    // unbond all parts
    for ( idx = range.first; idx < range.second; ++idx ) {
        x_BreakBond(idx);
    }
}


void CSeq_loc_CI_Impl::MakeBondA(size_t idx)
{
    pair<size_t, size_t> range;
    if ( IsInBond(idx) ) {
        range = GetBondRange(idx);
    }
    size_t bond_len = range.second-range.first;
    if ( bond_len > 0 ) {
        // update existing bond
        if ( idx != range.first ) {
            NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                           "CSeq_loc_I::MakeBondA(): "
                           "current position is B part of other bond");
        }
        if ( bond_len == 1 ) {
            // already bond with only part A
            return;
        }
        SetHasChanges();
        // unbond other parts
        for ( idx = range.first+1; idx < range.second; ++idx ) {
            x_BreakBond(idx);
        }
    }
    else {
        SetHasChanges();
        x_CreateBond(idx);
    }
}


void CSeq_loc_CI_Impl::MakeBondAB(size_t idx)
{
    if ( idx+1 >= m_Ranges.size() ) {
        NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                       "CSeq_loc_I::MakeBondAB(): "
                       "no more parts in the location");
    }
    pair<size_t, size_t> range;
    if ( IsInBond(idx) ) {
        range = GetBondRange(idx);
    }
    size_t bond_len = range.second-range.first;
    if ( bond_len > 0 ) {
        // update existing bond
        if ( idx != range.first ) {
            NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                           "CSeq_loc_I::MakeBondAB(): "
                           "current position is B part of other bond");
        }
        if ( bond_len == 2 ) {
            // already a bond with both parts A and B
            return;
        }
        SetHasChanges();
        if ( bond_len > 2 ) {
            // unbond extra parts after B
            for ( idx = range.first+2; idx < range.second; ++idx ) {
                x_BreakBond(idx);
            }
            return;
        }
        x_AddBondPartB(idx);
    }
    else {
        SetHasChanges();
        x_CreateBond(idx);
        x_AddBondPartB(idx);
    }
}


void CSeq_loc_CI_Impl::MakeBondB(size_t idx)
{
    if ( idx == 0 ) {
        NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                       "CSeq_loc_I::MakeBondB(): "
                       "no parts before current");
    }
    pair<size_t, size_t> range;
    if ( IsInBond(idx) ) {
        range = GetBondRange(idx);
    }
    else if ( IsInBond(idx-1) ) {
        range = GetBondRange(idx-1);
    }
    size_t bond_len = range.second-range.first;
    if ( bond_len > 0 ) {
        // update existing bond
        if ( range.first != idx + 1 ) {
            NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                           "CSeq_loc_I::MakeBondB(): "
                           "current position is not a B part of other bond");
        }
        if ( bond_len == 2 ) {
            // already a bond with both parts A and B
            return;
        }
        SetHasChanges();
        if ( bond_len > 2 ) {
            // unbond extra parts after B
            for ( idx = range.first+2; idx < range.second; ++idx ) {
                x_BreakBond(idx);
            }
            return;
        }
        _ASSERT(bond_len == 1);
        x_AddBondPartB(idx);
    }
    else {
        SetHasChanges();
        x_CreateBond(idx);
        x_AddBondPartB(idx);
    }
}


CRef<CSeq_loc> CSeq_loc_CI_Impl::MakeLocBond(size_t idx_begin,
                                             size_t idx_end) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    loc->SetBond().SetA(*MakePoint(m_Ranges[idx_begin]));
    if ( idx_begin+1 < idx_end ) {
        loc->SetBond().SetB(*MakePoint(m_Ranges[idx_begin+1]));
    }
    return loc;
}


static
CRef<CSeq_id> MakeId(const SSeq_loc_CI_RangeInfo& info)
{
    if ( !info.m_IdHandle ) {
        NCBI_THROW(CSeqLocException, eNotSet,
                   "CSeq_loc_I: part id is null");
    }
    return Ref(&const_cast<CSeq_id&>(*info.m_Id));
}


CRef<CInt_fuzz> CSeq_loc_CI_Impl::MakeFuzz(const CInt_fuzz& fuzz) const
{
    return Ref(&const_cast<CInt_fuzz&>(fuzz));
}


void CSeq_loc_CI_Impl::UpdatePoint(CSeq_point& pnt,
                                   const SSeq_loc_CI_RangeInfo& info) const
{
    pnt.SetId(*MakeId(info));
    pnt.SetPoint(info.m_Range.GetFrom());
    if ( info.m_IsSetStrand ) {
        pnt.SetStrand(info.m_Strand);
    }
    else {
        pnt.ResetStrand();
    }
    if ( info.m_Fuzz.first ) {
        pnt.SetFuzz(*MakeFuzz(*info.m_Fuzz.first));
    }
    else {
        pnt.ResetFuzz();
    }
}


void CSeq_loc_CI_Impl::SetPoint(SSeq_loc_CI_RangeInfo& info)
{
    info.m_Loc = MakeLocPoint(info);
}


CRef<CSeq_point>
CSeq_loc_CI_Impl::MakePoint(const SSeq_loc_CI_RangeInfo& info) const
{
    CRef<CSeq_point> seq_pnt(new CSeq_point);
    UpdatePoint(*seq_pnt, info);
    return seq_pnt;
}


CRef<CSeq_loc>
CSeq_loc_CI_Impl::MakeLocPoint(const SSeq_loc_CI_RangeInfo& info) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    loc->SetPnt(*MakePoint(info));
    return loc;
}


void CSeq_loc_CI_Impl::UpdatePoint(SSeq_loc_CI_RangeInfo& info)
{
    SetHasChanges();
    if ( WasPoint(info) ) {
        UpdatePoint(const_cast<CSeq_point&>(info.m_Loc->GetPnt()), info);
    }
}


void CSeq_loc_CI_Impl::UpdateLoc(SSeq_loc_CI_RangeInfo& info)
{
    SetHasChanges();
    if ( info.m_Loc ) {
        switch ( info.m_Loc->Which() ) {
        case CSeq_loc::e_Whole:
        case CSeq_loc::e_Empty:
        case CSeq_loc::e_Null:
        case CSeq_loc::e_Int:
        case CSeq_loc::e_Pnt:
            info.m_Loc = null;
            break;
        default:
            break;
        }
    }
}


CRef<CSeq_interval>
CSeq_loc_CI_Impl::MakeInterval(const SSeq_loc_CI_RangeInfo& info) const
{
    CRef<CSeq_interval> seq_int(new CSeq_interval);
    seq_int->SetId(*MakeId(info));
    seq_int->SetFrom(info.m_Range.GetFrom());
    seq_int->SetTo(info.m_Range.GetTo());
    if ( info.m_IsSetStrand ) {
        seq_int->SetStrand(info.m_Strand);
    }
    if ( info.m_Fuzz.first ) {
        seq_int->SetFuzz_from(*MakeFuzz(*info.m_Fuzz.first));
    }
    if ( info.m_Fuzz.second ) {
        seq_int->SetFuzz_to(*MakeFuzz(*info.m_Fuzz.second));
    }
    return seq_int;
}


CRef<CSeq_loc>
CSeq_loc_CI_Impl::MakeLocInterval(const SSeq_loc_CI_RangeInfo& info) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    loc->SetInt(*MakeInterval(info));
    return loc;
}


CRef<CSeq_loc>
CSeq_loc_CI_Impl::MakeLocOther(const SSeq_loc_CI_RangeInfo& info) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    if ( info.m_Range.IsWhole() ) {
        loc->SetWhole(*MakeId(info));
    }
    else {
        if ( !info.m_Range.Empty() ) {
            NCBI_THROW(CSeqLocException, eOtherError,
                       "CSeq_loc_I::MakeSeq_loc(): "
                       "cannot determine type of loc part");
        }
        if ( info.m_IdHandle ) {
            loc->SetEmpty(*MakeId(info));
        }
        else {
            loc->SetNull();
        }
    }
    return loc;
}


CRef<CSeq_loc>
CSeq_loc_CI_Impl::MakeRangeLoc(const SSeq_loc_CI_RangeInfo& info) const
{
    if ( ShouldBePoint(info) && CanBePoint(info) ) {
        return MakeLocPoint(info);
    }
    else if ( info.m_Range.IsWhole() || info.m_Range.Empty() ) {
        return MakeLocOther(info);
    }
    else {
        return MakeLocInterval(info);
    }
}


bool CSeq_loc_CI_Impl::CanBePacked_pnt(size_t idx_begin,
                                       size_t idx_end) const
{
    _ASSERT(idx_end >= idx_begin);
    if ( idx_end == idx_begin ) {
        return false;
    }
    const SSeq_loc_CI_RangeInfo& info0 = m_Ranges[idx_begin];
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        if ( IsInBond(idx) ) {
            return false;
        }
        const SSeq_loc_CI_RangeInfo& info = m_Ranges[idx];
        if ( !CanBePoint(m_Ranges[idx]) ) {
            return false;
        }
        if ( idx != idx_begin ) {
            if ( info.m_IdHandle != info0.m_IdHandle ) {
                return false;
            }
            if ( info.m_IsSetStrand != info0.m_IsSetStrand ) {
                return false;
            }
            if ( info0.m_IsSetStrand && info.m_Strand != info0.m_Strand ) {
                return false;
            }
            if ( info.m_Fuzz.first != info0.m_Fuzz.first ) {
                return false;
            }
        }
    }
    if ( HasEquivBreak(idx_begin, idx_end) ) {
        return false;
    }
    return true;
}


bool CSeq_loc_CI_Impl::CanBePacked_int(size_t idx_begin,
                                       size_t idx_end) const
{
    _ASSERT(idx_end >= idx_begin);
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        if ( IsInBond(idx) ) {
            return false;
        }
        if ( !CanBeInterval(idx) ) {
            return false;
        }
    }
    if ( HasEquivBreak(idx_begin, idx_end) ) {
        return false;
    }
    return true;
}


CRef<CSeq_loc> CSeq_loc_CI_Impl::MakeLocPacked_pnt(size_t idx_begin,
                                                   size_t idx_end) const
{
    const SSeq_loc_CI_RangeInfo& info0 = m_Ranges[idx_begin];
    CRef<CSeq_loc> loc(new CSeq_loc);
    CPacked_seqpnt& pnts = loc->SetPacked_pnt();
    pnts.SetId(*MakeId(info0));
    if ( info0.m_IsSetStrand ) {
        pnts.SetStrand(info0.m_Strand);
    }
    if ( info0.m_Fuzz.first ) {
        pnts.SetFuzz(*MakeFuzz(*info0.m_Fuzz.first));
    }
    pnts.SetPoints().reserve(idx_end-idx_begin);
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        const SSeq_loc_CI_RangeInfo& info = m_Ranges[idx];
        pnts.SetPoints().push_back(info.m_Range.GetFrom());
    }
    return loc;
}


CRef<CSeq_loc> CSeq_loc_CI_Impl::MakeLocPacked_int(size_t idx_begin,
                                                   size_t idx_end) const
{
    CRef<CSeq_loc> loc(new CSeq_loc);
    CPacked_seqint::Tdata& ints = loc->SetPacked_int().Set();
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        ints.push_back(MakeInterval(m_Ranges[idx]));
    }
    return loc;
}


static
void AddLoc(CRef<CSeq_loc>& loc, CRef<CSeq_loc> loc2)
{
    if ( !loc ) {
        loc = loc2;
        return;
    }
    if ( !loc->IsMix() ) {
        CRef<CSeq_loc> loc1 = loc;
        loc = new CSeq_loc;
        loc->SetMix().Set().push_back(loc1);
    }
    loc->SetMix().Set().push_back(loc2);
}


CRef<CSeq_loc> CSeq_loc_CI_Impl::MakeLoc(size_t idx_begin,
                                         size_t idx_end,
                                         CSeq_loc_I::EMakeType make_type,
                                         TUsedEquivs& used_equivs) const
{
    if ( const SEquivSet* eq_set = FindInnerEquivSet(idx_begin, idx_end, used_equivs) ) {
        used_equivs.insert(eq_set);
        size_t idx = eq_set->GetStartIndex();
        CRef<CSeq_loc> loc;
        if ( idx > idx_begin ) {
            loc = MakeLoc(idx_begin, idx, make_type, used_equivs);
        }
        CRef<CSeq_loc> loc2(new CSeq_loc);
        CSeq_loc_equiv& equiv = loc2->SetEquiv();
        ITERATE ( SEquivSet::TParts, it, eq_set->m_Parts ) {
            size_t part_end = eq_set->GetPartEndIndex(it);
            CRef<CSeq_loc> part_loc =
                MakeLoc(idx, part_end, make_type, used_equivs);
            equiv.Set().push_back(part_loc);
            idx = part_end;
        }
        AddLoc(loc, loc2);
        if ( idx < idx_end ) {
            loc2 = MakeLoc(idx, idx_end, make_type, used_equivs);
            if ( loc2->IsMix() ) {
                ITERATE ( CSeq_loc_mix::Tdata, it, loc2->GetMix().Get() ) {
                    AddLoc(loc, *it);
                }
            }
            else {
                AddLoc(loc, loc2);
            }
        }
        return loc;
    }

    if ( HasEquivBreak(idx_begin, idx_end) ) {
        NCBI_THROW(CSeqLocException, eBadLocation,
                   "CSeq_loc_I::MakeSeq_loc: equiv set overlaps with range");
    }

    if ( idx_begin == idx_end ) {
        CRef<CSeq_loc> loc(new CSeq_loc);
        loc->SetMix();
        return loc;
    }

    // check bonds
    const SSeq_loc_CI_RangeInfo& info0 = m_Ranges[idx_begin];
    if ( make_type == CSeq_loc_I::eMake_PreserveType && info0.m_Loc ) {
        // try to restore original CSeq_loc type
        if ( info0.m_Loc->IsPacked_pnt() &&
             CanBePacked_pnt(idx_begin, idx_end) ) {
            return MakeLocPacked_pnt(idx_begin, idx_end);
        }
        if ( info0.m_Loc->IsPacked_int() &&
             CanBePacked_int(idx_begin, idx_end) ) {
            return MakeLocPacked_int(idx_begin, idx_end);
        }
    }
    if ( make_type == CSeq_loc_I::eMake_CompactType &&
         idx_end - idx_begin > 1 ) {
        // try to make compact variants
        if ( CanBePacked_pnt(idx_begin, idx_end) ) {
            return MakeLocPacked_pnt(idx_begin, idx_end);
        }
        if ( CanBePacked_int(idx_begin, idx_end) ) {
            return MakeLocPacked_int(idx_begin, idx_end);
        }
    }

    CRef<CSeq_loc> loc;
    for ( size_t idx = idx_begin; idx < idx_end; ++idx ) {
        if ( IsInBond(idx) ) {
            pair<size_t, size_t> range = GetBondRange(idx);
            if ( range.first < idx || range.second > idx_end ||
                 !CanBeBond(range.first, range.second) ) {
                NCBI_THROW(CSeqLocException, eIncomatible,
                           "CSeq_loc_I::MakeSeq_loc: "
                           "invalid bond");
            }
            AddLoc(loc, MakeLocBond(range.first, range.second));
            idx = range.second-1;
            continue;
        }
        const SSeq_loc_CI_RangeInfo& info = m_Ranges[idx];
        if ( make_type == CSeq_loc_I::eMake_PreserveType ) {
            if ( WasPoint(info) && CanBePoint(info) ) {
                AddLoc(loc, MakeLocPoint(info));
                continue;
            }
            if ( (!info.m_Loc || info.m_Loc->IsInt()) &&
                 CanBeInterval(info) ) {
                AddLoc(loc, MakeLocInterval(info));
                continue;
            }
        }
        if ( CanBePoint(info) ) {
            AddLoc(loc, MakeLocPoint(info));
            continue;
        }
        if ( CanBeInterval(info) ) {
            AddLoc(loc, MakeLocInterval(info));
            continue;
        }
        AddLoc(loc, MakeLocOther(info));
    }
    if ( !loc ) {
        loc = new CSeq_loc;
    }
    return loc;
}


CRef<CSeq_loc> CSeq_loc_CI_Impl::MakeLoc(CSeq_loc_I::EMakeType make_type) const
{
    TUsedEquivs used_equivs;
    return MakeLoc(0, m_Ranges.size(), make_type, used_equivs);
}


CSeq_loc_CI::CSeq_loc_CI(void)
    : m_Impl(new CSeq_loc_CI_Impl),
      m_Index(0)
{
}


CSeq_loc_CI::CSeq_loc_CI(const CSeq_loc& loc,
                         EEmptyFlag empty_flag,
                         ESeqLocOrder order)
    : m_Impl(new CSeq_loc_CI_Impl(loc, empty_flag, order)),
      m_Index(0)
{
}


CSeq_loc_CI::CSeq_loc_CI(const CSeq_loc_CI& iter, size_t pos)
    : m_Impl(iter.m_Impl),
      m_Index(0)
{
    SetPos(pos);
}


CSeq_loc_CI::~CSeq_loc_CI()
{
}


CSeq_loc_CI::CSeq_loc_CI(const CSeq_loc_CI& iter)
    : m_Impl(iter.m_Impl),
      m_Index(iter.m_Index)
{
}


CSeq_loc_CI& CSeq_loc_CI::operator= (const CSeq_loc_CI& iter)
{
    m_Impl = iter.m_Impl;
    m_Index = iter.m_Index;
    return *this;
}


bool CSeq_loc_CI::operator== (const CSeq_loc_CI& iter) const
{
    // Check if both are at end.
    bool is_end = m_Impl->IsEnd(m_Index);
    if ( iter.m_Impl->IsEnd(iter.m_Index) ) {
        return is_end;
    }
    else {
        return !is_end && m_Impl == iter.m_Impl && m_Index == iter.m_Index;
    }
}


bool CSeq_loc_CI::operator!= (const CSeq_loc_CI& iter) const
{
    return !(*this == iter);
}


bool CSeq_loc_CI::IsInBond(void) const
{
    x_CheckValid("IsInBond()");
    return m_Impl->IsInBond(m_Index);
}


bool CSeq_loc_CI::IsBondA(void) const
{
    x_CheckValid("IsBondA()");
    return m_Impl->IsBondPartA(m_Index);
}


bool CSeq_loc_CI::IsBondB(void) const
{
    x_CheckValid("IsBondB()");
    return m_Impl->IsBondPartB(m_Index);
}


pair<CSeq_loc_CI, CSeq_loc_CI>
CSeq_loc_CI::GetBondRange(void) const
{
    x_CheckValid("GetBondRange()");
    pair<size_t, size_t> indexes = m_Impl->GetBondRange(m_Index);
    return make_pair(CSeq_loc_CI(*this, indexes.first),
                     CSeq_loc_CI(*this, indexes.second));
}


bool CSeq_loc_CI::HasEquivSets(void) const
{
    return m_Impl->HasEquivSets();
}


bool CSeq_loc_CI::IsInEquivSet(void) const
{
    x_CheckValid("IsInEquivSet()");
    return m_Impl->IsInEquivSet(m_Index);
}


size_t CSeq_loc_CI::GetEquivSetsCount(void) const
{
    x_CheckValid("GetEquivSetsCount()");
    return m_Impl->GetEquivSetsCount(m_Index);
}


pair<CSeq_loc_CI, CSeq_loc_CI>
CSeq_loc_CI::GetEquivSetRange(size_t level) const
{
    x_CheckValid("GetEquivSetRange()");
    pair<size_t, size_t> indexes = m_Impl->GetEquivSetRange(m_Index, level);
    return make_pair(CSeq_loc_CI(*this, indexes.first),
                     CSeq_loc_CI(*this, indexes.second));
}


pair<CSeq_loc_CI, CSeq_loc_CI>
CSeq_loc_CI::GetEquivPartRange(size_t level) const
{
    x_CheckValid("GetEquivPartRange()");
    pair<size_t, size_t> indexes = m_Impl->GetEquivPartRange(m_Index, level);
    return make_pair(CSeq_loc_CI(*this, indexes.first),
                     CSeq_loc_CI(*this, indexes.second));
}


const CSeq_loc& CSeq_loc_CI::GetSeq_loc(void) const
{
    return GetEmbeddingSeq_loc();
}


const CSeq_loc& CSeq_loc_CI::GetEmbeddingSeq_loc(void) const
{
    x_CheckValid("GetEmbeddingSeq_loc()");
    CConstRef<CSeq_loc> loc = x_GetRangeInfo().m_Loc;
    if ( !loc ) {
        NCBI_THROW(CSeqLocException, eNotSet,
                   "CSeq_loc_CI::GetSeq_loc(): NULL seq-loc");
    }
    return *loc;
}


CConstRef<CSeq_loc> CSeq_loc_CI::GetRangeAsSeq_loc(void) const
{
    x_CheckValid("GetRangeAsSeq_loc()");
    const SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Loc ) {
        switch ( info.m_Loc->Which() ) {
            // Single-range, empty or whole seq-loc can be used as-is.
        case CSeq_loc::e_not_set:
        case CSeq_loc::e_Null:
        case CSeq_loc::e_Empty:
        case CSeq_loc::e_Whole:
        case CSeq_loc::e_Int:
        case CSeq_loc::e_Pnt:
            return info.m_Loc;
        default:
            break;
        }
    }
    // Create a new seq-loc for the current range.
    CRef<CSeq_loc> loc = m_Impl->MakeRangeLoc(info);
    return ConstRef(loc.Release());
}


bool CSeq_loc_CI::x_IsValid(void) const
{
    return m_Impl && m_Index < m_Impl->GetRanges().size();
}


void CSeq_loc_CI::x_ThrowNotValid(const char* where) const
{
    NCBI_THROW_FMT(CSeqLocException, eBadIterator,
                   x_GetIteratorType() << "::" << where << ": "
                   "iterator is not valid");
}


const char* CSeq_loc_CI::x_GetIteratorType(void) const
{
    return "CSeq_loc_CI";
}


const SSeq_loc_CI_RangeInfo& CSeq_loc_CI::x_GetRangeInfo(void) const
{
    // The index validity must be checked by the caller.
    return m_Impl->GetRanges()[m_Index];
}


size_t CSeq_loc_CI::GetSize(void) const
{
    return m_Impl->GetRanges().size();
}


void CSeq_loc_CI::SetPos(size_t pos)
{
    if ( pos > GetSize() ) {
        NCBI_THROW_FMT(CSeqLocException, eOtherError,
                       x_GetIteratorType() << "::SetPos(): "
                       "position is too big: "<<pos<<" > "<<GetSize());
    }
    m_Index = pos;
}


/////////////////////////////////////////////////////////////////////////////
// CSeq_loc_I

CSeq_loc_I::CSeq_loc_I(void)
{
}


CSeq_loc_I::CSeq_loc_I(CSeq_loc& loc)
    : CSeq_loc_CI(loc, eEmpty_Allow, eOrder_Biological)
{
}


CSeq_loc_I::CSeq_loc_I(const CSeq_loc_I& iter, size_t pos)
    : CSeq_loc_CI(iter)
{
    SetPos(pos);
}


CSeq_loc_I::~CSeq_loc_I()
{
}


bool CSeq_loc_I::x_IsValidForInsert(void) const
{
    return m_Impl && m_Index <= m_Impl->GetRanges().size();
}


const char* CSeq_loc_I::x_GetIteratorType(void) const
{
    return "CSeq_loc_I";
}


SSeq_loc_CI_RangeInfo& CSeq_loc_I::x_GetRangeInfo(void)
{
    // The index validity must be checked by the caller.
    return m_Impl->GetRanges()[m_Index];
}


void CSeq_loc_I::x_SetSeq_id_Handle(SSeq_loc_CI_RangeInfo& info,
                                    const CSeq_id_Handle& id)
{
    info.m_Id = id.GetSeqId();
    info.m_IdHandle = id;
}


bool CSeq_loc_I::HasChanges(void) const
{
    return m_Impl->HasChanges();
}


CSeq_loc_I::EEquivMode CSeq_loc_I::GetEquivMode(void) const
{
    return m_Impl->GetEquivMode();
}


void CSeq_loc_I::SetEquivMode(EEquivMode mode)
{
    m_Impl->SetEquivMode(mode);
}


void CSeq_loc_I::Delete(void)
{
    x_CheckValid("Delete()");
    m_Impl->DeleteRange(m_Index);
}


CSeq_loc_I CSeq_loc_I::InsertNull(void)
{
    x_CheckValidForInsert("InsertNull()");
    m_Impl->InsertRange(m_Index, CSeq_loc::e_Null);
    return CSeq_loc_I(*this, m_Index++);
}


CSeq_loc_I CSeq_loc_I::InsertEmpty(const CSeq_id_Handle& id)
{
    x_CheckValidForInsert("InsertEmpty()");
    SSeq_loc_CI_RangeInfo& info =
        m_Impl->InsertRange(m_Index, CSeq_loc::e_Empty);
    x_SetSeq_id_Handle(info, id);
    return CSeq_loc_I(*this, m_Index++);
}


CSeq_loc_I CSeq_loc_I::InsertWhole(const CSeq_id_Handle& id)
{
    x_CheckValidForInsert("InsertWhole()");
    SSeq_loc_CI_RangeInfo& info =
        m_Impl->InsertRange(m_Index, CSeq_loc::e_Whole);
    x_SetSeq_id_Handle(info, id);
    info.m_Range = TRange::GetWhole();
    return CSeq_loc_I(*this, m_Index++);
}


CSeq_loc_I CSeq_loc_I::InsertPoint(const CSeq_id_Handle& id,
                                   TSeqPos pos,
                                   ENa_strand strand)
{
    x_CheckValidForInsert("InsertPoint()");
    SSeq_loc_CI_RangeInfo& info =
        m_Impl->InsertRange(m_Index, CSeq_loc::e_Pnt);
    x_SetSeq_id_Handle(info, id);
    info.m_Range.SetFrom(pos).SetLength(1);
    if ( strand != eNa_strand_unknown ) {
        info.SetStrand(strand);
    }
    m_Impl->SetPoint(info);
    return CSeq_loc_I(*this, m_Index++);
}


CSeq_loc_I CSeq_loc_I::InsertInterval(const CSeq_id_Handle& id,
                                      const TRange& range,
                                      ENa_strand strand)
{
    x_CheckValidForInsert("InsertInterval()");
    SSeq_loc_CI_RangeInfo& info =
        m_Impl->InsertRange(m_Index, CSeq_loc::e_Int);
    x_SetSeq_id_Handle(info, id);
    info.m_Range = range;
    if ( strand != eNa_strand_unknown ) {
        info.SetStrand(strand);
    }
    if ( !range.IsWhole() && !range.Empty() ) {
        info.m_Loc = m_Impl->MakeLocInterval(info);
    }
    return CSeq_loc_I(*this, m_Index++);
}


void CSeq_loc_I::SetSeq_id_Handle(const CSeq_id_Handle& id)
{
    x_CheckValid("SetSeq_id_Handle()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_IdHandle != id ) {
        x_SetSeq_id_Handle(info, id);
        m_Impl->UpdatePoint(info);
    }
}


void CSeq_loc_I::SetRange(const TRange& range)
{
    x_CheckValid("SetRange()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Range != range ) {
        info.m_Range = range;
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::SetFrom(TSeqPos from)
{
    x_CheckValid("SetFrom()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Range.GetFrom() != from ) {
        info.m_Range.SetFrom(from);
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::SetTo(TSeqPos to)
{
    x_CheckValid("SetTo()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Range.GetTo() != to ) {
        info.m_Range.SetTo(to);
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::SetPoint(TSeqPos pos)
{
    x_CheckValid("SetPoint()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( !m_Impl->WasPoint(info) || info.m_Range != TRange(pos, pos) ) {
        info.m_Range = TRange(pos, pos);
        if ( m_Impl->WasPoint(info) ) {
            m_Impl->UpdatePoint(info);
        }
        else {
            m_Impl->SetPoint(info);
        }
    }
}


void CSeq_loc_I::ResetStrand(void)
{
    x_CheckValid("ResetStrand()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_IsSetStrand ) {
        info.m_IsSetStrand = false;
        info.m_Strand = eNa_strand_unknown;
        m_Impl->UpdatePoint(info);
    }
}


void CSeq_loc_I::SetStrand(ENa_strand strand)
{
    x_CheckValid("SetStrand()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( !info.m_IsSetStrand || info.m_Strand != strand ) {
        info.SetStrand(strand);
        m_Impl->UpdatePoint(info);
    }
}


void CSeq_loc_I::ResetFuzzFrom(void)
{
    x_CheckValid("ResetFuzzFrom()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Fuzz.first ) {
        info.m_Fuzz.first = null;
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::SetFuzzFrom(CInt_fuzz& fuzz)
{
    x_CheckValid("SetFuzzFrom()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( !info.m_Fuzz.first || !info.m_Fuzz.first->Equals(fuzz) ) {
        info.m_Fuzz.first = SerialClone(fuzz);
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::ResetFuzzTo(void)
{
    x_CheckValid("ResetFuzzTo()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Fuzz.second ) {
        info.m_Fuzz.second = null;
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::SetFuzzTo(CInt_fuzz& fuzz)
{
    x_CheckValid("SetFuzzTo()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( !info.m_Fuzz.second || !info.m_Fuzz.second->Equals(fuzz) ) {
        info.m_Fuzz.second = SerialClone(fuzz);
        m_Impl->UpdateLoc(info);
    }
}


void CSeq_loc_I::ResetFuzz(void)
{
    x_CheckValid("ResetFuzz()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( info.m_Fuzz.first || info.m_Fuzz.second ) {
        info.m_Fuzz.first = info.m_Fuzz.second = null;
        m_Impl->UpdatePoint(info);
    }
}


void CSeq_loc_I::SetFuzz(CInt_fuzz& fuzz)
{
    x_CheckValid("SetFuzz()");
    SSeq_loc_CI_RangeInfo& info = x_GetRangeInfo();
    if ( !info.m_Fuzz.first || !info.m_Fuzz.first->Equals(fuzz) ||
         info.m_Fuzz.first != info.m_Fuzz.second ) {
        info.m_Fuzz.first = info.m_Fuzz.second = SerialClone(fuzz);
        m_Impl->UpdatePoint(info);
    }
}


CRef<CSeq_loc> CSeq_loc_I::MakeSeq_loc(EMakeType make_type) const
{
    return m_Impl->MakeLoc(make_type);
}


pair<CSeq_loc_I, CSeq_loc_I>
CSeq_loc_I::GetEquivSetRange(size_t level) const
{
    x_CheckValid("GetEquivSetRange()");
    pair<size_t, size_t> indexes = m_Impl->GetEquivSetRange(m_Index, level);
    return make_pair(CSeq_loc_I(*this, indexes.first),
                     CSeq_loc_I(*this, indexes.second));
}


pair<CSeq_loc_I, CSeq_loc_I>
CSeq_loc_I::GetEquivPartRange(size_t level) const
{
    x_CheckValid("GetEquivPartRange()");
    pair<size_t, size_t> indexes = m_Impl->GetEquivPartRange(m_Index, level);
    return make_pair(CSeq_loc_I(*this, indexes.first),
                     CSeq_loc_I(*this, indexes.second));
}


void CSeq_loc_I::RemoveBond(void)
{
    x_CheckValid("RemoveBond()");
    m_Impl->RemoveBond(m_Index);
}


void CSeq_loc_I::MakeBondA(void)
{
    x_CheckValid("MakeBondA()");
    m_Impl->MakeBondA(m_Index);
}


void CSeq_loc_I::MakeBondAB(void)
{
    x_CheckValid("MakeBondAB()");
    m_Impl->MakeBondAB(m_Index);
}


void CSeq_loc_I::MakeBondB(void)
{
    x_CheckValid("MakeBondB()");
    m_Impl->MakeBondB(m_Index);
}


// CSeq_loc_I
/////////////////////////////////////////////////////////////////////////////


// Append a string representation of a CSeq_id to label
inline
void s_GetLabel(const CSeq_id& id, string* label)
{
    CNcbiOstrstream os;
    id.WriteAsFasta(os);
    *label += CNcbiOstrstreamToString(os);
}


// Append to label info for a CSeq_point
inline
const CSeq_id* s_GetLabel
(const CSeq_point& pnt,
 const CSeq_id*    last_id,
 string*           label)
{
    // If CSeq_id does not match last_id, then append id to label
    if ( !last_id  ||  !last_id->Match(pnt.GetId()) ) {
        s_GetLabel(pnt.GetId(), label);
        *label += ":";
    }

    // Add strand info to label
    if (pnt.IsSetStrand()) {
        *label += GetTypeInfo_enum_ENa_strand()
            ->FindName(pnt.GetStrand(), true);
    }

    if (pnt.IsSetFuzz()) {
        // Append Fuzz info to label
        pnt.GetFuzz().GetLabel(label, pnt.GetPoint());
    } else {
        // Append 1 based point to label
        *label += NStr::IntToString(pnt.GetPoint()+1);
    }

    // update last_id
    last_id = &pnt.GetId();
    return last_id;
}


// Append to label info for CSeq_interval
inline
const CSeq_id* s_GetLabel
(const CSeq_interval& itval,
 const CSeq_id*       last_id,
 string*              label)
{
    if (!last_id || !last_id->Match(itval.GetId())) {
        s_GetLabel(itval.GetId(), label);
        *label += ":";
    }
    last_id = &itval.GetId();
    //if (itval.IsSetStrand() && itval.GetStrand() != eNa_strand_unknown) {
    if (itval.IsSetStrand() && (itval.GetStrand() == eNa_strand_minus || itval.GetStrand() == eNa_strand_both_rev)) {
        //*label += GetTypeInfo_enum_ENa_strand()->FindName(itval.GetStrand(), true);
        *label += "c";
    }
    if (itval.IsSetStrand() &&
        (itval.GetStrand() == eNa_strand_minus ||
         itval.GetStrand() == eNa_strand_both_rev)) {
        if (itval.IsSetFuzz_to()) {
            itval.GetFuzz_to().GetLabel(label, itval.GetTo(), false);
        } else {
            *label += NStr::IntToString(itval.GetTo()+1);
        }
        *label += "-";
        if (itval.IsSetFuzz_from()) {
            itval.GetFuzz_from().GetLabel(label, itval.GetFrom());
        } else {
            *label += NStr::IntToString(itval.GetFrom()+1);
        }
    } else {
        if (itval.IsSetFuzz_from()) {
            itval.GetFuzz_from().GetLabel
                (label, itval.GetFrom(), false);
        } else {
            *label += NStr::IntToString(itval.GetFrom()+1);
        }
        *label += "-";
        if (itval.IsSetFuzz_to()) {
            itval.GetFuzz_to().GetLabel(label, itval.GetTo());
        } else {
            *label += NStr::IntToString(itval.GetTo()+1);
        }
    }
    return last_id;
}


// Forward declaration
const CSeq_id* s_GetLabel
(const CSeq_loc& loc,
 const CSeq_id*  last_id,
 string*         label,
 bool            first = false);


// Appends to label info for each CSeq_loc in a list of CSeq_locs
inline
const CSeq_id* s_GetLabel
(const list<CRef<CSeq_loc> >&  loc_list,
 const CSeq_id*                last_id,
 string*                       label)
{
    bool first = true;
    ITERATE (list<CRef<CSeq_loc> >, it, loc_list) {

        // Append to label for each CSeq_loc in list
        last_id = s_GetLabel(**it, last_id, label, first);
        first = false;
    }

    return last_id;
}


// Builds a label based upon a CSeq_loc and all embedded CSeq_locs
const CSeq_id* s_GetLabel
(const CSeq_loc& loc,
 const CSeq_id*  last_id,
 string*         label,
 bool            first)
{
    // Ensure label is not null
    if (!label) {
        return last_id;
    }

    // Put a comma separator if necessary
    if (!first) {
        *label += ", ";
    }

    // Loop through embedded CSeq_locs and create a label, depending on
    // type of CSeq_loc
    switch (loc.Which()) {
    case CSeq_loc::e_Null:
        *label += "~";
        break;
    case CSeq_loc::e_Empty:
        *label += "{";
        s_GetLabel(loc.GetEmpty(), label);
        last_id = &loc.GetEmpty();
        *label += "}";
        break;
    case CSeq_loc::e_Whole:
        s_GetLabel(loc.GetWhole(), label);
        last_id = &loc.GetWhole();
        break;
    case CSeq_loc::e_Int:
        last_id = s_GetLabel(loc.GetInt(), last_id, label);
        break;
    case CSeq_loc::e_Packed_int:
    {
        *label += "(";
        bool frst = true;
        ITERATE(CPacked_seqint::Tdata, it, loc.GetPacked_int().Get()) {
            if (!frst) {
                *label += ", ";
            }
            frst = false;
            last_id = s_GetLabel(**it, last_id, label);
        }
        *label += ")";
        break;
    }
    case CSeq_loc::e_Pnt:
        last_id = s_GetLabel(loc.GetPnt(), last_id, label);
        break;
    case CSeq_loc::e_Packed_pnt:
        *label += "(" + loc.GetPacked_pnt().GetId().AsFastaString() + ":";
        {{
             string str;
             ITERATE (CPacked_seqpnt::TPoints, iter,
                      loc.GetPacked_pnt().GetPoints()) {
                 if ( !str.empty() ) {
                     str += ", ";
                 }
                 str += NStr::IntToString(*iter);
             }
             *label += str;
         }}
        *label += ")";
        last_id = &loc.GetPacked_pnt().GetId();
        break;
    case CSeq_loc::e_Mix:
        *label += "[";
        last_id = s_GetLabel(loc.GetMix().Get(), last_id, label);
        *label += "]";
        break;
    case CSeq_loc::e_Equiv:
        *label += "[";
        last_id = s_GetLabel(loc.GetEquiv().Get(), last_id, label);
        *label += "]";
        break;
    case CSeq_loc::e_Bond:
        last_id = s_GetLabel(loc.GetBond().GetA(), last_id, label);
        *label += "=";
        if (loc.GetBond().IsSetB()) {
            last_id = s_GetLabel(loc.GetBond().GetB(), last_id, label);
        } else {
            *label += "?";
        }
        break;
    case CSeq_loc::e_Feat:
        *label += "(feat)";
        break;
    default:
        *label += "(?\?)";
        break;
    }
    return last_id;
}


bool CSeq_loc::IsPartialStart(ESeqLocExtremes ext) const
{
    switch (Which ()) {
        case e_Null :
            return true;

        case e_Int :
            return GetInt().IsPartialStart(ext);

        case e_Packed_int :
            return GetPacked_int().IsPartialStart(ext);

        case e_Pnt :
            return GetPnt().IsPartialStart(ext);

        case e_Packed_pnt :
            return GetPacked_pnt().IsPartialStart(ext);

        case e_Mix :
            return GetMix().IsPartialStart(ext);

        default :
            break;
    }

    return false;
}


bool CSeq_loc::IsPartialStop(ESeqLocExtremes ext) const
{
    switch (Which ()) {
        case e_Null :
            return true;

        case e_Int :
            return GetInt().IsPartialStop(ext);

        case e_Packed_int :
            return GetPacked_int().IsPartialStop(ext);

        case e_Pnt :
            return GetPnt().IsPartialStop(ext);

        case e_Packed_pnt :
            return GetPacked_pnt().IsPartialStop(ext);

        case e_Mix :
            return GetMix().IsPartialStop(ext);

        default :
            break;
    }

    return false;
}


void CSeq_loc::SetPartialStart(bool val, ESeqLocExtremes ext)
{
    if (val == IsPartialStart(ext)) {
        return;
    }

    switch (Which()) {
        case e_Int:
            SetInt().SetPartialStart(val, ext);
            break;

        case e_Packed_int :
            SetPacked_int().SetPartialStart(val, ext);
            break;

        case e_Pnt:
            SetPnt().SetPartialStart(val, ext);
            break;

        case e_Packed_pnt:
            SetPacked_pnt().SetPartialStart(val, ext);
            break;

        case e_Mix :
            SetMix().SetPartialStart(val, ext);
            break;

        default :
            break;
    }
}


void CSeq_loc::SetPartialStop(bool val, ESeqLocExtremes ext)
{
    if (val == IsPartialStop(ext)) {
        return;
    }

    switch (Which()) {
        case e_Int:
            SetInt().SetPartialStop(val, ext);
            break;

        case e_Packed_int :
            SetPacked_int().SetPartialStop(val, ext);
            break;

        case e_Pnt:
            SetPnt().SetPartialStop(val, ext);
            break;

        case e_Packed_pnt:
            SetPacked_pnt().SetPartialStop(val, ext);
            break;

        case e_Mix:
            SetMix().SetPartialStop(val, ext);
            break;

        default :
            break;
    }
}


bool CSeq_loc::IsTruncatedStart(ESeqLocExtremes ext) const
{
    switch (Which ()) {
        case e_Int :
            return GetInt().IsTruncatedStart(ext);

        case e_Packed_int :
            return GetPacked_int().IsTruncatedStart(ext);

        case e_Pnt :
            return GetPnt().IsTruncatedStart(ext);

        case e_Packed_pnt :
            return GetPacked_pnt().IsTruncatedStart(ext);

        case e_Mix :
            return GetMix().IsTruncatedStart(ext);

        default :
            break;
    }

    return false;
}


bool CSeq_loc::IsTruncatedStop(ESeqLocExtremes ext) const
{
    switch (Which ()) {
        case e_Int :
            return GetInt().IsTruncatedStop(ext);

        case e_Packed_int :
            return GetPacked_int().IsTruncatedStop(ext);

        case e_Pnt :
            return GetPnt().IsTruncatedStop(ext);

        case e_Packed_pnt :
            return GetPacked_pnt().IsTruncatedStop(ext);

        case e_Mix :
            return GetMix().IsTruncatedStop(ext);

        default :
            break;
    }

    return false;
}


void CSeq_loc::SetTruncatedStart(bool val, ESeqLocExtremes ext)
{
    if (val == IsTruncatedStart(ext)) {
        return;
    }

    switch (Which()) {
        case e_Int:
            SetInt().SetTruncatedStart(val, ext);
            break;

        case e_Packed_int :
            SetPacked_int().SetTruncatedStart(val, ext);
            break;

        case e_Pnt:
            SetPnt().SetTruncatedStart(val, ext);
            break;

        case e_Packed_pnt:
            SetPacked_pnt().SetTruncatedStart(val, ext);
            break;

        case e_Mix :
            SetMix().SetTruncatedStart(val, ext);
            break;

        default :
            break;
    }
}


void CSeq_loc::SetTruncatedStop(bool val, ESeqLocExtremes ext)
{
    if (val == IsTruncatedStop(ext)) {
        return;
    }

    switch (Which()) {
        case e_Int:
            SetInt().SetTruncatedStop(val, ext);
            break;

        case e_Packed_int :
            SetPacked_int().SetTruncatedStop(val, ext);
            break;

        case e_Pnt:
            SetPnt().SetTruncatedStop(val, ext);
            break;

        case e_Packed_pnt:
            SetPacked_pnt().SetTruncatedStop(val, ext);
            break;

        case e_Mix:
            SetMix().SetTruncatedStop(val, ext);
            break;

        default :
            break;
    }
}


// Appends a label suitable for display (e.g., error messages)
// label must point to an existing string object
// Method just returns if label is null
void CSeq_loc::GetLabel(string* label) const
{
    s_GetLabel(*this, 0, label, true);
}


// assign the 'id' field of each sub-interval to the supplied id
void CSeq_loc::SetId(CSeq_id& id)
{
    InvalidateCache();
    switch (Which()) {
    case e_Null:
        break;

    case e_Int:
        SetInt().SetId(id);
        break;

    case e_Pnt:
        SetPnt().SetId(id);
        break;

    case e_Packed_int:
        NON_CONST_ITERATE (CPacked_seqint::Tdata, iter, SetPacked_int().Set()) {
            (*iter)->SetId(id);
        }
        break;

    case e_Packed_pnt:
        SetPacked_pnt().SetId(id);
        break;

    case e_Mix:
        NON_CONST_ITERATE (CSeq_loc_mix::Tdata, iter, SetMix().Set()) {
            (*iter)->SetId(id);
        }
        break;

    case e_Whole:
        SetWhole(id);
        break;

    case e_Empty:
        SetEmpty(id);
        break;

    case e_Equiv:
        NON_CONST_ITERATE (CSeq_loc_equiv::Tdata, iter, SetEquiv().Set()) {
            (*iter)->SetId(id);
        }
        break;

    case e_Bond:
        if (GetBond().IsSetA()) {
            SetBond().SetA().SetId(id);
        }
        if (GetBond().IsSetB()) {
            SetBond().SetB().SetId(id);
        }
        break;

    case e_Feat:
        ERR_POST_X(1, Error
                      << "unhandled loc type in CSeq_loc::SetId(): e_Feat");
        break;

    default:
        ERR_POST_X(2, Error << "unhandled loc type in CSeq_loc::SetId(): "
                      << Which());
        break;
    }
}


bool CSeq_loc::x_CheckId(const CSeq_id*& id, bool may_throw) const
{
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
        {
            // no Seq-id
            return true;
        }
    case e_Empty:
        {
            return x_UpdateId(id, &GetEmpty(), may_throw);
        }
    case e_Whole:
        {
            return x_UpdateId(id, &GetWhole(), may_throw);
        }
    case e_Int:
        {
            const CSeq_interval& loc = GetInt();
            return x_UpdateId(id, &loc.GetId(), may_throw);
        }
    case e_Pnt:
        {
            const CSeq_point& pnt = GetPnt();
            return x_UpdateId(id, &pnt.GetId(), may_throw);
        }
    case e_Packed_int:
        {
            // Check ID of each interval
            const CPacked_seqint& ints = GetPacked_int();
            ITERATE ( CPacked_seqint::Tdata, ii, ints.Get() ) {
                const CSeq_interval& loc = **ii;
                if ( !x_UpdateId(id, &loc.GetId(), may_throw) ) {
                    return false;
                }
            }
            return true;
        }
    case e_Packed_pnt:
        {
            const CPacked_seqpnt& pnts = GetPacked_pnt();
            return x_UpdateId(id, &pnts.GetId(), may_throw);
            break;
        }
    case e_Mix:
        {
            // Check ID of each sub-location.
            const CSeq_loc_mix& mix = GetMix();
            ITERATE( CSeq_loc_mix::Tdata, li, mix.Get() ) {
                if ( !(*li)->CheckId(id, may_throw) ) {
                    return false;
                }
            }
            return true;
        }
    case e_Bond:
        {
            const CSeq_bond& bond = GetBond();
            if ( bond.CanGetA() ) {
                if ( !x_UpdateId(id, &bond.GetA().GetId(), may_throw) ) {
                    return false;
                }
            }
            if ( bond.CanGetB() ) {
                return x_UpdateId(id, &bond.GetB().GetId(), may_throw);
            }
            return true;
        }
    case e_Equiv:
        {
            // Doesn't make much sense to test equiv, but ...
            ITERATE(CSeq_loc_equiv::Tdata, li, GetEquiv().Get()) {
                if ( !(*li)->CheckId(id, may_throw) ) {
                    return false;
                }
            }
            return true;
        }
    case e_Feat:
    default:
        if (may_throw) {
            NCBI_THROW_FMT(CSeqLocException, eUnsupported,
                           "CSeq_loc::CheckId(): "
                           "unsupported location type: "<<
                           SelectionName(Which()));
        }
        return false;
    }
}


void CSeq_loc::ChangeToMix(void)
{
    switch ( Which() ) {
    case e_not_set:
        {
            SetMix();
            break;
        }
    case e_Mix:
        {
            break;
        }
    case e_Packed_int:
        {
            // unpack
            CRef<CSeq_loc> self(new CSeq_loc);
            self->Assign(*this, eShallow);

            CSeq_loc_mix& mix = SetMix();
            NON_CONST_ITERATE (CPacked_seqint::Tdata, it, self->SetPacked_int().Set()) {
                CRef<CSeq_loc> ival(new CSeq_loc);
                ival->SetInt(**it);
                mix.Set().push_back(ival);
            }
            break;
        }
    default:
        {
            CRef<CSeq_loc> self(new CSeq_loc);
            self->Assign(*this, eShallow);
            CSeq_loc_mix& mix = SetMix();
            mix.AddSeqLoc(*self);
        }
    }
}


void CSeq_loc::ChangeToPackedInt(void)
{
    switch ( Which() ) {
    case e_not_set:
    case e_Null:
        {
            SetPacked_int();
            return;
        }
    case e_Packed_int:
        {
            return;
        }
    case e_Int:
        {
            CConstRef<CSeq_interval> self(&GetInt());
            SetPacked_int().AddInterval(*self);
            return;
        }
    case e_Pnt:
        {
            // Make an interval of length 1
            CRef<CSeq_interval> new_int(new CSeq_interval);
            new_int->SetId().Assign(GetPnt().GetId());
            new_int->SetFrom(GetPnt().GetPoint());
            new_int->SetTo(GetPnt().GetPoint());
            if (GetPnt().IsSetStrand()) {
                new_int->SetStrand(GetPnt().GetStrand());
            }
            if (GetPnt().IsSetFuzz()) {
                const CInt_fuzz& fuzz = GetPnt().GetFuzz();
                if (!fuzz.IsLim()  ||  fuzz.GetLim() != CInt_fuzz::eLim_gt) {
                    new_int->SetFuzz_from().Assign(fuzz);
                }
                if (!fuzz.IsLim()  ||  fuzz.GetLim() != CInt_fuzz::eLim_lt) {
                    new_int->SetFuzz_to().Assign(fuzz);
                }
            }
            SetPacked_int().AddInterval(*new_int);
            return;
        }
    case e_Mix:
        {
            // Recursively convert each sub-location to packed-int, then merge.
            // Work with copies of sub-locs so that if an exception is thrown
            // the original location remains unchanged.
            vector<CRef<CSeq_loc> > sub_locs;
            sub_locs.reserve(GetMix().Get().size());
            ITERATE (CSeq_loc_mix::Tdata, orig_sub_loc, GetMix().Get()) {
                CRef<CSeq_loc> new_sub_loc(new CSeq_loc);
                new_sub_loc->Assign(**orig_sub_loc);
                new_sub_loc->ChangeToPackedInt();
                sub_locs.push_back(new_sub_loc);
            }
            SetPacked_int();  // in case there are zero intervals
            ITERATE (vector<CRef<CSeq_loc> >, sub_loc, sub_locs) {
                copy((*sub_loc)->GetPacked_int().Get().begin(),
                     (*sub_loc)->GetPacked_int().Get().end(),
                     back_inserter(SetPacked_int().Set()));
            }
            return;
        }
    default:
        NCBI_THROW_FMT(CSeqLocException, eIncomatible,
                       "CSeq_loc::ChangeToPackedInt(): "
                       "Can not convert location to packed-int: "<<
                       SelectionName(Which()));
    }
}


void CSeq_loc::x_ChangeToMix(const CSeq_loc& other)
{
    ChangeToMix();
    SetMix().AddSeqLoc(other);
}


void CSeq_loc::x_ChangeToPackedInt(const CSeq_interval& other)
{
    _ASSERT(IsInt());

    ChangeToPackedInt();
    SetPacked_int().AddInterval(other);
}


void CSeq_loc::x_ChangeToPackedInt(const CSeq_loc& other)
{
    _ASSERT(IsInt());
    _ASSERT(other.IsInt()  ||  other.IsPacked_int());

    ChangeToPackedInt();

    if ( other.IsInt() ) {
        SetPacked_int().AddInterval(other.GetInt());
    } else {  // other is packed int
        SetPacked_int().AddIntervals(other.GetPacked_int());
    }
}


void CSeq_loc::x_ChangeToPackedPnt(const CSeq_loc& other)
{
    _ASSERT(IsPnt());
    _ASSERT(other.IsPnt()  ||  other.IsPacked_pnt());

    CRef<CSeq_point> pnt(&SetPnt());
    CPacked_seqpnt& ppnt = SetPacked_pnt();
    if ( pnt->IsSetStrand() ) {
        ppnt.SetStrand(pnt->GetStrand());
    }
    if ( pnt->IsSetId() ) {
        ppnt.SetId(pnt->SetId());
    }
    if ( pnt->IsSetFuzz() ) {
        ppnt.SetFuzz(pnt->SetFuzz());
    }
    ppnt.AddPoint(pnt->GetPoint());

    if ( other.IsPnt() ) {
        ppnt.AddPoint(other.GetPnt().GetPoint());
    } else { // other is packed point
        ppnt.AddPoints(other.GetPacked_pnt().GetPoints());
    }
}


template<typename T1, typename T2>
bool s_CanAdd(const T1& obj1, const T2& obj2)
{
    // test strands
    {{
        ENa_strand s1 = obj1.CanGetStrand() ? obj1.GetStrand() : eNa_strand_unknown;
        ENa_strand s2 = obj2.CanGetStrand() ? obj2.GetStrand() : eNa_strand_unknown;
        if ( s1 != s2 ) {
            return false;
        }
    }}

    // test ids
    {{
        const CSeq_id* id1 = obj1.CanGetId() ? &obj1.GetId() : 0;
        const CSeq_id* id2 = obj2.CanGetId() ? &obj2.GetId() : 0;
        // both null - ok to add (this will probably never happen)
        if (id1 == 0  &&  id2 == 0) return true;
        // ids are different (one may be null)
        if (id1 == 0  ||  id2 == 0  ||  !id1->Match(*id2) ) {
            return false;
        }
    }}

    // test fuzz
    {{
        const CInt_fuzz* f1 = obj1.CanGetFuzz() ? &obj1.GetFuzz() : 0;
        const CInt_fuzz* f2 = obj2.CanGetFuzz() ? &obj2.GetFuzz() : 0;
        // both null - ok to add
        if (f1 == 0  &&  f2 == 0) return true;
        // fuzzes are different (one may be null)
        if (f1 == 0  ||  f2 == 0  ||  !f1->Equals(*f2)) {
            return false;
        }
    }}

    return true;
}


bool s_CanAdd(const CSeq_loc& loc1, const CSeq_loc& loc2)
{
    switch ( loc1.Which() ) {
    case CSeq_loc::e_Pnt:
        {
            switch ( loc2.Which() ) {
            case CSeq_loc::e_Pnt:
                return s_CanAdd(loc1.GetPnt(), loc2.GetPnt());
            case CSeq_loc::e_Packed_pnt:
                return s_CanAdd(loc1.GetPnt(), loc2.GetPacked_pnt());
            default:
                break;
            }
            break;
        }
    case CSeq_loc::e_Packed_pnt:
        {
            switch ( loc2.Which() ) {
            case CSeq_loc::e_Pnt:
                return s_CanAdd(loc1.GetPacked_pnt(), loc2.GetPnt());
            case CSeq_loc::e_Packed_pnt:
                return s_CanAdd(loc1.GetPacked_pnt(), loc2.GetPacked_pnt());
            default:
                break;
            }
            break;
        }
    default:
        {
            return false;
        }
    }

    return false;
}


void CSeq_loc::Add(const CSeq_loc& other)
{
    InvalidateCache();
    switch ( Which() ) {
    case e_not_set:
        {
            Assign(other);
            break;
        }
    case e_Null:
        {
            // ??? skip if other is null?
            x_ChangeToMix(other);
            break;
        }
    case e_Empty:
        {
            // ??? skip if other is empty and ids match?
            x_ChangeToMix(other);
            break;
        }

    case e_Whole:
        {
            x_ChangeToMix(other);
            break;
        }
    case e_Int:
        {
            if ( other.IsInt()  ||  other.IsPacked_int() ) {
                x_ChangeToPackedInt(other);
            } else {
                x_ChangeToMix(other);
            }
        }
        break;
    case e_Pnt:
        {
            if ( s_CanAdd(*this, other) ) {
                x_ChangeToPackedPnt(other);
            } else {
                x_ChangeToMix(other);
            }
            break;
        }
    case e_Packed_int:
        {
            if ( other.IsInt() ) {
                SetPacked_int().AddInterval(other.GetInt());
            } else if ( other.IsPacked_int() ) {
                SetPacked_int().AddIntervals(other.GetPacked_int());
            } else {
                x_ChangeToMix(other);
            }
            break;
        }
    case e_Packed_pnt:
        {
            if ( s_CanAdd(*this, other) ) {
                if ( other.IsPnt() ) {
                    SetPacked_pnt().AddPoint(other.GetPnt().GetPoint());
                } else if ( other.IsPacked_pnt() ) {
                    SetPacked_pnt().AddPoints(other.GetPacked_pnt().GetPoints());
                }
            } else {
                x_ChangeToMix(other);
            }
            break;
        }
    case e_Mix:
        {
            SetMix().AddSeqLoc(other);
            break;
        }
    case e_Equiv:
        {
            SetEquiv().Add(other);
            break;
        }
    case e_Bond:
        {
            x_ChangeToMix(other);
            break;
        }
    case e_Feat:
    default:
        NCBI_THROW_FMT(CSeqLocException, eIncomatible,
                       "CSeq_loc::Add(): "
                       "cannot add sub-location to location of type: "<<
                       SelectionName(Which()));
    }
}


void CSeq_loc::FlipStrand(void)
{
    switch ( Which() ) {
    case e_Int:
        SetInt().FlipStrand();
        break;
    case e_Pnt:
        SetPnt().FlipStrand();
        break;
    case e_Packed_int:
        SetPacked_int().FlipStrand();
        break;
    case e_Packed_pnt:
        SetPacked_pnt().FlipStrand();
        break;
    case e_Mix:
        SetMix().FlipStrand();
        break;

    default:
        break;
    }
}


// Types used in operations with seq-locs

class CRangeWithFuzz : public CSeq_loc::TRange
{
public:
    typedef CSeq_loc::TRange TParent;
    typedef CConstRef<CInt_fuzz>  TFuzz;

    CRangeWithFuzz(const TParent& rg)
        : TParent(rg), m_Strand(eNa_strand_unknown)
    {
    }
    CRangeWithFuzz(const CSeq_loc_CI& it)
        : TParent(it.GetRange()),
          m_Fuzz_from(it.GetFuzzFrom()),
          m_Fuzz_to(it.GetFuzzTo()),
          m_Strand(it.GetStrand())
    {
    }

    void ResetFuzzFrom(void) { m_Fuzz_from.Reset(); }
    void ResetFuzzTo(void) { m_Fuzz_to.Reset(); }
    bool IsSetFuzzFrom(void) const { return m_Fuzz_from; }
    bool IsSetFuzzTo(void) const { return m_Fuzz_to; }
    const CInt_fuzz& GetFuzzFrom(void) const { return *m_Fuzz_from; }
    const CInt_fuzz& GetFuzzTo(void) const { return *m_Fuzz_to; }

    // Add fuzzes assuming that both ranges had the same 'from'
    void AddFuzzFrom(const CRangeWithFuzz& rg)
    {
        x_AddFuzz(m_Fuzz_from, rg.m_Fuzz_from, rg.m_Strand);
    }

    // Add fuzzes assuming that both ranges had the same 'to'
    void AddFuzzTo(const CRangeWithFuzz& rg)
    {
        x_AddFuzz(m_Fuzz_to, rg.m_Fuzz_to, rg.m_Strand);
    }

    void AddFuzzFrom(const CInt_fuzz& fuzz, ENa_strand strand)
    {
        x_AddFuzz(m_Fuzz_from, ConstRef(&fuzz), strand);
    }
    
    void AddFuzzTo(const CInt_fuzz& fuzz, ENa_strand strand)
    {
        x_AddFuzz(m_Fuzz_to, ConstRef(&fuzz), strand);
    }

    void CopyFrom(const CRangeWithFuzz& rg)
    {
        SetFrom(rg.GetFrom());
        m_Fuzz_from = rg.m_Fuzz_from;
    }

    void CopyTo(const CRangeWithFuzz& rg)
    {
        SetTo(rg.GetTo());
        m_Fuzz_to = rg.m_Fuzz_to;
    }

    CRangeWithFuzz& operator +=(const CRangeWithFuzz& rg)
    {
        TParent::position_type old_from = GetFrom();
        TParent::position_type old_to = GetTo();
        TParent::operator+=(rg);
        if (old_from != GetFrom()) {
            m_Fuzz_from.Reset(rg.m_Fuzz_from);
        }
        else if (old_from == rg.GetFrom()) {
            // Reset fuzz if it's not the same for both ranges
            AddFuzzFrom(rg);
        }
        if (old_to != GetTo()) {
            m_Fuzz_to.Reset(rg.m_Fuzz_to);
        }
        else if (old_to == rg.GetTo()) {
            AddFuzzTo(rg);
        }
        return *this;
    }

    CRangeWithFuzz& operator +=(const TParent& rg)
    {
        TParent::position_type old_from = GetFrom();
        TParent::position_type old_to = GetTo();
        TParent::operator+=(rg);
        // Reset fuzz if the corresponding extreme changes
        if (old_from != GetFrom()) {
            ResetFuzzFrom();
        }
        if (old_to != GetTo()) {
            ResetFuzzTo();
        }
        return *this;
    }

private:
    CRef<CInt_fuzz> x_SetFuzz(TFuzz&           fuzz,
                              const CInt_fuzz* copy_from)
    {
        TFuzz copy_from_cref;
        if (copy_from == fuzz) copy_from_cref.Reset(copy_from);
        // Since TFuzz is a const-ref, setting fuzz requires creating
        // a new object
        CRef<CInt_fuzz> new_fuzz(new CInt_fuzz);
        // The new value is optional
        if ( copy_from ) {
            new_fuzz->Assign(*copy_from);
        }
        fuzz.Reset(new_fuzz);
        return new_fuzz;
    }

    void x_AddFuzz(TFuzz&       fuzz,
                   const TFuzz& other,
                   ENa_strand   other_strand)
    {
        if ( !fuzz ) {
            // Use fuzz from the other range if available
            if ( other ) {
                x_SetFuzz(fuzz, other.GetPointerOrNull());
            }
            return;
        }
        if ( !other ) {
            // The other range has no fuzz, keep the current one
            return;
        }
        if (fuzz->Which() != other->Which()) {
            // Fuzzes have different types, reset to lim-unk.
            CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, NULL);
            new_fuzz->SetLim(CInt_fuzz::eLim_unk);
            return;
        }

        const CInt_fuzz& fz = *fuzz;
        const CInt_fuzz& ofz = *other;
        // Both fuzzes are set and have the same type, try to merge them
        switch ( fz.Which() ) {
        case CInt_fuzz::e_Lim:
            {
                CInt_fuzz::ELim this_lim = fz.GetLim();
                CInt_fuzz::ELim other_lim = ofz.GetLim();
                bool this_rev = IsReverse(m_Strand);
                bool other_rev = IsReverse(other_strand);
                bool other_lt = other_lim == CInt_fuzz::eLim_lt  ||
                    (!other_rev  &&  other_lim == CInt_fuzz::eLim_tl)  ||
                    (other_rev  &&  other_lim == CInt_fuzz::eLim_tr);
                bool other_gt = other_lim == CInt_fuzz::eLim_gt  ||
                    (!other_rev  &&  other_lim == CInt_fuzz::eLim_tr)  ||
                    (other_rev  &&  other_lim == CInt_fuzz::eLim_tl);
                switch ( fz.GetLim() ) {
                case CInt_fuzz::eLim_lt:
                    if ( other_lt ) {
                        return; // the same
                    }
                    break;
                case CInt_fuzz::eLim_gt:
                    if ( other_gt ) {
                        return; // the same
                    }
                    break;
                case CInt_fuzz::eLim_tl:
                    if ((!this_rev  &&  other_lt)  ||
                        (this_rev  &&  other_gt)) {
                        return; // the same
                    }
                    break;
                case CInt_fuzz::eLim_tr:
                    if ((!this_rev  &&  other_gt)  ||
                        (this_rev  &&  other_lt)) {
                        return; // the same
                    }
                    break;
                default:
                    if (other_lim == this_lim) {
                        return;
                    }
                    break;
                }
                // Different limits - reset to lim-unk.
                CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, NULL);
                new_fuzz->SetLim(CInt_fuzz::eLim_unk);
                break;
            }
        case CInt_fuzz::e_Alt:
            {
                // Use union
                CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, fuzz);
                new_fuzz->SetAlt().insert(
                    new_fuzz->SetAlt().end(),
                    ofz.GetAlt().begin(),
                    ofz.GetAlt().end());
                break;
            }
        case CInt_fuzz::e_Range:
            {
                // Use union
                CInt_fuzz::C_Range::TMin min1 = fz.GetRange().GetMin();
                CInt_fuzz::C_Range::TMin min2 = ofz.GetRange().GetMin();
                CInt_fuzz::C_Range::TMax max1 = fz.GetRange().GetMax();
                CInt_fuzz::C_Range::TMax max2 = ofz.GetRange().GetMax();
                if (min1 > min2  ||  max1 < max2) {
                    CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, NULL);
                    new_fuzz->SetRange().SetMin(min1 < min2 ? min1 : min2);
                    new_fuzz->SetRange().SetMax(max1 > max2 ? max1 : max2);
                }
                break;
            }
        case CInt_fuzz::e_P_m:
            {
                // Use max value
                CInt_fuzz::TP_m pm = ofz.GetP_m();
                if (fz.GetP_m() < pm) {
                    CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, NULL);
                    new_fuzz->SetP_m(pm);
                }
                break;
            }
        case CInt_fuzz::e_Pct:
            {
                // Use max value
                CInt_fuzz::TPct pct = ofz.GetPct();
                if (fz.GetPct() < pct) {
                    CRef<CInt_fuzz> new_fuzz = x_SetFuzz(fuzz, NULL);
                    new_fuzz->SetPct(pct);
                }
                break;
            }
        default:
            // Failed to merge fuzzes
            fuzz.Reset();
            break;
        }
    }

    TFuzz m_Fuzz_from;
    TFuzz m_Fuzz_to;
    ENa_strand m_Strand;
};


class CSeq_id_Handle_Wrapper
{
public:
    CSeq_id_Handle_Wrapper(void) {}

    CSeq_id_Handle_Wrapper(const CSeq_id_Handle& idh, const CSeq_id& id)
        : m_Handle(idh)
    {
        if (id.IsLocal() && id.GetLocal().IsStr()) {
            m_Id.Reset(&id);
        }
    }

    CConstRef<CSeq_id> GetSeqId(void) const { return m_Id ? m_Id : m_Handle.GetSeqId(); }

    const CSeq_id_Handle& GetHandle(void) const { return m_Handle; }

    bool operator== (const CSeq_id_Handle_Wrapper& handle) const
        {
            return m_Handle == handle.m_Handle;
        }

    bool operator!= (const CSeq_id_Handle_Wrapper& handle) const
        {
            return m_Handle != handle.m_Handle;
        }
    bool operator<  (const CSeq_id_Handle_Wrapper& handle) const
        {
            return m_Handle < handle.m_Handle;
        }

    DECLARE_OPERATOR_BOOL(m_Handle);

private:
    CSeq_id_Handle m_Handle;
    CConstRef<CSeq_id> m_Id;
};


typedef CRangeWithFuzz                          TRangeWithFuzz;
typedef vector<TRangeWithFuzz>                  TRanges;
typedef map<CSeq_id_Handle_Wrapper, TRanges>    TIdToRangeMap;
typedef CRangeCollection<TSeqPos>               TRangeColl;
typedef map<CSeq_id_Handle_Wrapper, TRangeColl> TIdToRangeColl;


class CRange_Less
{
public:
    CRange_Less(void) {}

    bool operator() (const TRangeWithFuzz& rg1, const TRangeWithFuzz& rg2) const
    {
        if ( rg1.IsWhole() ) {
            return !rg2.IsWhole();
        }
        if ( rg1.Empty() ) {
            return !rg2.Empty()  &&  !rg2.IsWhole();
        }
        return !rg2.IsWhole()  &&  !rg2.Empty()  &&  rg1 < rg2;
    }
};


class CRange_ReverseLess
{
public:
    CRange_ReverseLess(void) {}

    bool operator() (const TRangeWithFuzz& rg1, const TRangeWithFuzz& rg2) const
    {
        if ( rg1.IsWhole() ) {
            return !rg2.IsWhole();
        }
        if ( rg1.Empty() ) {
            return !rg2.Empty()  &&  !rg2.IsWhole();
        }
        if ( rg2.IsWhole()  ||  rg2.Empty() ) {
            return false;
        }
        if (rg1.GetTo() != rg2.GetTo()) {
            return rg1.GetTo() > rg2.GetTo();
        }
        return rg1.GetFrom() < rg2.GetFrom();
    }
};


static inline
bool x_MatchStrand(ENa_strand str1, ENa_strand str2, CSeq_loc::TOpFlags flags)
{
    return (flags & CSeq_loc::fStrand_Ignore) != 0 ||
    IsReverse(str1) == IsReverse(str2);
}


static
bool x_MergeRanges(TRangeWithFuzz& rg1, ENa_strand str1,
                   const TRangeWithFuzz& rg2, ENa_strand str2,
                   CSeq_loc::TOpFlags flags)
{
    if ( !x_MatchStrand(str1, str2, flags) ) {
        return false;
    }
    // Check contained
    if ( (flags & CSeq_loc::fMerge_Contained) != 0 ) {
        if (rg1.GetFrom() <= rg2.GetFrom()  &&  rg1.GetTo() >= rg2.GetTo()) {
            // rg2 already contained in rg1
            if (rg1.GetFrom() == rg2.GetFrom()) {
                rg1.AddFuzzFrom(rg2);
            }
            if (rg1.GetTo() == rg2.GetTo()) {
                rg1.AddFuzzTo(rg2);
            }
            return true;
        }
        if (rg1.GetFrom() >= rg2.GetFrom()  &&  rg1.GetTo() <= rg2.GetTo()) {
            // rg1 contained in rg2
            bool same_from = rg1.GetFrom() == rg2.GetFrom();
            bool same_to = rg1.GetTo() == rg2.GetTo();
            rg1 = rg2;
            if (same_from) {
                rg1.AddFuzzFrom(rg2);
            }
            if (same_to) {
                rg1.AddFuzzTo(rg2);
            }
            return true;
        }
    }
    // Check overlapping
    if ( (flags & CSeq_loc::fMerge_OverlappingOnly) != 0  &&
        rg1.IntersectingWith(rg2) ) {
        rg1 += rg2;
        return true;
    }
    // Check abutting
    if ((flags & CSeq_loc::fMerge_AbuttingOnly) != 0) {
        if ( !IsReverse(str1) ) {
            if ( rg1.GetToOpen() == rg2.GetFrom() ) {
                rg1.CopyTo(rg2);
                return true;
            }
        }
        else {
            if (rg1.GetFrom() == rg2.GetToOpen()) {
                rg1.CopyFrom(rg2);
                return true;
            }
        }
    }
    return false;
}


static
void x_PushRange(CSeq_loc& dst,
                 const CSeq_id_Handle_Wrapper& idh,
                 const TRangeWithFuzz& rg,
                 ENa_strand strand)
{
    if (dst.Which() != CSeq_loc::e_not_set) {
        if ( !dst.IsMix() ) {
            dst.ChangeToMix();
        }
    }
    if ( !idh ) {
        // NULL
        if (dst.IsMix()) {
            dst.SetMix().Set().push_back(Ref(new CSeq_loc(CSeq_loc::e_Null)));
        }
        else {
            dst.SetNull();
        }
        return;
    }
    CRef<CSeq_id> id(new CSeq_id());
    id->Assign(*idh.GetSeqId());
    if ( rg.IsWhole() ) {
        if (dst.IsMix()) {
            CRef<CSeq_loc> whole(new CSeq_loc);
            whole->SetWhole(*id);
            dst.SetMix().Set().push_back(whole);
        }
        else {
            dst.SetWhole(*id);
        }
    }
    else if ( rg.Empty() ) {
        if (dst.IsMix()) {
            CRef<CSeq_loc> empty(new CSeq_loc);
            empty->SetEmpty(*id);
            dst.SetMix().Set().push_back(empty);
        }
        else {
            dst.SetEmpty(*id);
        }
    }
    else if ( rg.GetLength() == 1 &&
        rg.IsSetFuzzFrom() == rg.IsSetFuzzTo() &&
        ( !rg.IsSetFuzzFrom() ||
            rg.GetFuzzFrom().Equals(rg.GetFuzzTo()) ) )
    {
        // Preserve points
        CRef<CSeq_point> pnt(new CSeq_point);
        pnt->SetId(*id);
        pnt->SetPoint(rg.GetFrom());
        if (strand != eNa_strand_unknown) {
            pnt->SetStrand(strand);
        }
        if ( rg.IsSetFuzzFrom() ) {
            pnt->SetFuzz().Assign(rg.GetFuzzFrom());
        } else if( rg.IsSetFuzzTo() ) {
            pnt->SetFuzz().Assign(rg.GetFuzzTo());
        }
        if (dst.IsMix()) {
            CRef<CSeq_loc> pnt_loc(new CSeq_loc);
            pnt_loc->SetPnt(*pnt);
            dst.SetMix().Set().push_back(pnt_loc);
        }
        else {
            dst.SetPnt(*pnt);
        }
    }
    else {
        if (dst.IsMix()) {
            CRef<CSeq_loc> int_loc(new CSeq_loc);
            CSeq_interval& ival = int_loc->SetInt();
            ival.SetFrom(rg.GetFrom());
            ival.SetTo(rg.GetTo());
            ival.SetId().Assign(*id);
            if (strand != eNa_strand_unknown) {
                ival.SetStrand(strand);
            }
            if ( rg.IsSetFuzzFrom() ) {
                ival.SetFuzz_from().Assign(rg.GetFuzzFrom());
            }
            if ( rg.IsSetFuzzTo() ) {
                ival.SetFuzz_to().Assign(rg.GetFuzzTo());
            }
            dst.SetMix().Set().push_back(int_loc);
        }
        else {
            CRef<CSeq_interval> interval(new CSeq_interval(*id,
                rg.GetFrom(),
                rg.GetTo(),
                strand));
            if ( rg.IsSetFuzzFrom() ) {
                interval->SetFuzz_from().Assign(rg.GetFuzzFrom());
            }
            if ( rg.IsSetFuzzTo() ) {
                interval->SetFuzz_to().Assign(rg.GetFuzzTo());
            }
            dst.SetInt(*interval);
        }
    }
}


static
void x_SingleRange(CSeq_loc& dst,
                   const CSeq_loc& src,
                   ISynonymMapper& syn_mapper)
{
    // Create a single range
    TRangeWithFuzz total_rg(TRangeWithFuzz::GetEmpty());
    CSeq_id_Handle_Wrapper first_id;
    ENa_strand first_strand = eNa_strand_unknown;
    for (CSeq_loc_CI it(src, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper next_id(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        if ( !next_id ) {
            // Ignore NULLs
            continue;
        }
        if ( first_id ) {
            // Seq-id may be missing for NULL seq-loc
            if (next_id  &&  first_id != next_id) {
                NCBI_THROW(CSeqLocException, eMultipleId,
                           "Can not merge multi-id seq-loc");
            }
        }
        else {
            first_id = next_id;
            first_strand = it.GetStrand();
        }
        total_rg += TRangeWithFuzz(it);
    }
    if ( first_id ) {
        CRef<CSeq_id> id(new CSeq_id);
        id->Assign(*first_id.GetSeqId());
        CRef<CSeq_interval> interval(new CSeq_interval(*id,
                                                       total_rg.GetFrom(),
                                                       total_rg.GetTo(),
                                                       first_strand));
        if ( total_rg.IsSetFuzzFrom() ) {
            interval->SetFuzz_from().Assign(total_rg.GetFuzzFrom());
        }
        if ( total_rg.IsSetFuzzTo() ) {
            interval->SetFuzz_to().Assign(total_rg.GetFuzzTo());
        }
        dst.SetInt(*interval);
    }
    else {
        // Null seq-loc
        dst.SetNull();
    }
}


static
void x_RangesToSeq_loc(CSeq_loc& dst,
                       TIdToRangeMap& id_map,
                       ENa_strand default_strand,
                       CSeq_loc::TOpFlags flags)
{
    // Iterate ids for each strand
    NON_CONST_ITERATE(TIdToRangeMap, id_it, id_map) {
        if ( !id_it->first ) {
            // All NULLs merged
            x_PushRange(dst,
                        id_it->first,
                        TRangeWithFuzz(TRangeWithFuzz::GetEmpty()),
                        eNa_strand_unknown);
            continue;
        }
        CRef<CSeq_id> id(new CSeq_id);
        id->Assign(*id_it->first.GetSeqId());
        TRanges& ranges = id_it->second;
        if ( (flags & CSeq_loc::fSort) != 0 ) {
            if ( !IsReverse(default_strand) ) {
                sort(ranges.begin(), ranges.end(), CRange_Less());
            }
            else {
                sort(ranges.begin(), ranges.end(), CRange_ReverseLess());
            }
        }
        // Merge ranges according to the flags, add to destination
        TRangeWithFuzz last_rg(TRangeWithFuzz::GetEmpty());
        bool have_range = false;
        ITERATE(TRanges, rg, ranges) {
            if (x_MergeRanges(last_rg, default_strand,
                            *rg, default_strand,
                            flags)) {
                have_range = true;
                continue;
            }
            // No merging - push current range, reset last values
            if (have_range) {
                x_PushRange(dst, id_it->first, last_rg, default_strand);
            }
            last_rg = *rg;
            have_range = true;
        }
        if (have_range) {
            x_PushRange(dst, id_it->first, last_rg, default_strand);
        }
    }
}


static
void x_MergeNoSort(CSeq_loc& dst,
                   const CSeq_loc& src,
                   CSeq_loc::TOpFlags flags,
                   ISynonymMapper& syn_mapper)
{
    _ASSERT((flags & CSeq_loc::fSort) == 0);
    CSeq_id_Handle_Wrapper last_id;
    TRangeWithFuzz last_rg(TRangeWithFuzz::GetEmpty());
    ENa_strand last_strand = eNa_strand_unknown;
    bool have_range = false;
    for (CSeq_loc_CI it(src, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper idh(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        // ID and strand must match
        TRangeWithFuzz it_rg(it);
        if ( have_range  &&  last_id == idh ) {
            if (x_MergeRanges(last_rg,
                              last_strand,
                              it_rg,
                              it.GetStrand(),
                              flags)) {
                have_range = true;
                continue;
            }
        }
        // No merging - push current range, reset last values
        if ( have_range ) {
            x_PushRange(dst, last_id, last_rg, last_strand);
        }
        last_id = idh;
        last_rg = it_rg;
        last_strand = it.GetStrand();
        have_range = true;
    }
    if ( have_range ) {
        x_PushRange(dst, last_id, last_rg, last_strand);
    }
    if (dst.Which() == CSeq_loc::e_not_set) {
        dst.SetNull();
    }
}


static
void x_MergeAndSort(CSeq_loc& dst,
                    const CSeq_loc& src,
                    CSeq_loc::TOpFlags flags,
                    ISynonymMapper& syn_mapper)
{
    bool use_strand = (flags & CSeq_loc::fStrand_Ignore) == 0;
    // Id -> range map for both strands
    unique_ptr<TIdToRangeMap> pid_map_minus(use_strand ?
        new TIdToRangeMap : 0);
    TIdToRangeMap id_map_plus;
    TIdToRangeMap& id_map_minus = use_strand ?
        *pid_map_minus.get() : id_map_plus;

    // Prepare default strands
    ENa_strand default_plus = eNa_strand_unknown;
    ENa_strand default_minus = eNa_strand_unknown;

    // Split location by by id/strand/range
    for (CSeq_loc_CI it(src, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper idh(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        ENa_strand strand = it.GetStrand();
        if ( IsReverse(strand) ) {
            id_map_minus[idh].push_back(TRangeWithFuzz(it));
            if (use_strand) {
                if (default_minus == eNa_strand_unknown) {
                    default_minus = strand;
                }
                else if (default_minus != strand) {
                    default_minus = eNa_strand_minus;
                }
            }
        }
        else {
            id_map_plus[idh].push_back(TRangeWithFuzz(it));
            if (use_strand) {
                if (default_plus == eNa_strand_unknown) {
                    default_plus = strand;
                }
                else if (default_plus != strand) {
                    default_plus = eNa_strand_plus;
                }
            }
        }
    }

    x_RangesToSeq_loc(dst, id_map_plus, default_plus, flags);
    if ( use_strand ) {
        x_RangesToSeq_loc(dst, id_map_minus, default_minus, flags);
    }
    if (dst.Which() == CSeq_loc::e_not_set) {
        dst.SetNull();
    }
}


typedef map<TSeqPos, CConstRef<CInt_fuzz> > TFuzzMap;

static
void x_SingleRange(CSeq_loc& dst,
                   const CSeq_loc& minuend,
                   const TFuzzMap& fuzz_from_plus,
                   const TFuzzMap& fuzz_from_minus,
                   const TFuzzMap& fuzz_to_plus,
                   const TFuzzMap& fuzz_to_minus,
                   TIdToRangeColl& rg_coll_plus,
                   TIdToRangeColl& rg_coll_minus,
                   ISynonymMapper& syn_mapper,
                   ILengthGetter& len_getter)
{
    TRangeWithFuzz total_rg(TRangeWithFuzz::GetEmpty());
    CSeq_id_Handle_Wrapper first_id;
    ENa_strand first_strand = eNa_strand_unknown;
    for (CSeq_loc_CI it(minuend, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper next_id(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        if ( !next_id ) {
            // Ignore NULLs
            continue;
        }
        if ( first_id ) {
            // Seq-id may be missing for NULL seq-loc
            if (next_id  &&  first_id != next_id) {
                NCBI_THROW(CSeqLocException, eMultipleId,
                           "Can not merge multi-id seq-loc");
            }
        }
        else {
            first_id = next_id;
            first_strand = it.GetStrand();
        }
        TRangeWithFuzz it_range = TRangeWithFuzz(it);
        if (it_range.GetFrom() >= total_rg.GetFrom()  &&
            it_range.GetTo() <= total_rg.GetTo()) {
            // Nothing new can be added from this interval
            continue;
        }
        if ( it_range.IsWhole() ) {
            it_range.SetOpen(0, len_getter.GetLength(it.GetSeq_id()));
            it_range.ResetFuzzFrom();
            it_range.ResetFuzzTo();
        }
        TRangeColl it_rg_coll(it_range);
        TIdToRangeColl& rg_coll = IsReverse(it.GetStrand()) ?
            rg_coll_minus : rg_coll_plus;
        TIdToRangeColl::const_iterator id_it = rg_coll.find(next_id);
        if (id_it != rg_coll.end()) {
            it_rg_coll -= id_it->second;
        }
        TRangeWithFuzz curr_rg(it_rg_coll.GetLimits());
        if (curr_rg.GetFrom() == it_range.GetFrom()) {
            curr_rg.AddFuzzFrom(it_range);
        }
        if (curr_rg.GetTo() == it_range.GetTo()) {
            curr_rg.AddFuzzTo(it_range);
        }
        
        // If range start/stop comes from subtrahend, copy fuzz.
        const TFuzzMap& fm_from = IsReverse(it.GetStrand()) ? fuzz_from_minus : fuzz_from_plus;
        TFuzzMap::const_iterator subtr_fuzz = fm_from.find(curr_rg.GetToOpen());
        if (subtr_fuzz != fm_from.end()) {
            curr_rg.AddFuzzTo(*subtr_fuzz->second, it.GetStrand());
        }
        const TFuzzMap& fm_to = IsReverse(it.GetStrand()) ? fuzz_to_minus : fuzz_to_plus;
        subtr_fuzz = fm_to.find(curr_rg.GetFrom());
        if (subtr_fuzz != fm_to.end()) {
            curr_rg.AddFuzzFrom(*subtr_fuzz->second, it.GetStrand());
        }
            
        total_rg += curr_rg;
    }

    if ( first_id ) {
        CRef<CSeq_id> id(new CSeq_id);
        id->Assign(*first_id.GetSeqId());
        CRef<CSeq_interval> interval(new CSeq_interval(*id,
                                                       total_rg.GetFrom(),
                                                       total_rg.GetTo(),
                                                       first_strand));
        if ( total_rg.IsSetFuzzFrom() ) {
            CRef<CInt_fuzz> fuzz(new CInt_fuzz);
            fuzz->Assign(total_rg.GetFuzzFrom());
            interval->SetFuzz_from(*fuzz);
        }
        if ( total_rg.IsSetFuzzTo() ) {
            CRef<CInt_fuzz> fuzz(new CInt_fuzz);
            fuzz->Assign(total_rg.GetFuzzTo());
            interval->SetFuzz_to(*fuzz);
        }
        dst.SetInt(*interval);
    }
    else {
        // Null seq-loc
        dst.SetNull();
    }
}


static
void x_SubNoSort(CSeq_loc& dst,
                 const CSeq_loc& minuend,
                 const TFuzzMap& fuzz_from_plus,
                 const TFuzzMap& fuzz_from_minus,
                 const TFuzzMap& fuzz_to_plus,
                 const TFuzzMap& fuzz_to_minus,
                 TIdToRangeColl& rg_coll_plus,
                 TIdToRangeColl& rg_coll_minus,
                 ISynonymMapper& syn_mapper,
                 ILengthGetter& len_getter,
                 CSeq_loc::TOpFlags flags)
{
    _ASSERT((flags & CSeq_loc::fSort) == 0);
    CSeq_id_Handle_Wrapper last_id;
    TRangeWithFuzz last_rg(TRangeWithFuzz::GetEmpty());
    ENa_strand last_strand = eNa_strand_unknown;
    bool have_range = false;
    for (CSeq_loc_CI it(minuend, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper idh(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        bool rev = IsReverse(it.GetStrand());
        TRangeWithFuzz it_range = TRangeWithFuzz(it);
        if ( it_range.IsWhole() ) {
            it_range.SetOpen(0, len_getter.GetLength(it.GetSeq_id()));
        }
        TRangeColl it_rg_coll(it_range);
        TIdToRangeColl& rg_coll = IsReverse(it.GetStrand()) ?
            rg_coll_minus : rg_coll_plus;
        TIdToRangeColl::const_iterator id_it = rg_coll.find(idh);
        list<TRangeWithFuzz> result_ranges;
        bool modified = false;
        if (id_it != rg_coll.end()) {
            // Check if there's anything to subtract
            ITERATE(TRangeColl, check_it, id_it->second) {
                if ( it_rg_coll.IntersectingWith(*check_it) ) {
                    it_rg_coll -= id_it->second;
                    modified = true;
                    ITERATE(TRangeColl, result_it, it_rg_coll) {
                        if ( rev ) {
                            result_ranges.push_front(*result_it);
                        }
                        else {
                            result_ranges.push_back(*result_it);
                        }
                    }
                    break;
                }
            }
        }
        if ( modified ) {
            ITERATE(list<TRangeWithFuzz>, rg_it, result_ranges) {
                TRangeWithFuzz curr_rg(*rg_it);
                if (curr_rg.GetFrom() == it_range.GetFrom()) {
                    curr_rg.AddFuzzFrom(it_range);
                }
                if (curr_rg.GetTo() == it_range.GetTo()) {
                    curr_rg.AddFuzzTo(it_range);
                }

                // If range start/stop comes from subtrahend, copy fuzz.
                const TFuzzMap& fm_from = IsReverse(it.GetStrand()) ? fuzz_from_minus : fuzz_from_plus;
                TFuzzMap::const_iterator subtr_fuzz = fm_from.find(curr_rg.GetToOpen());
                if (subtr_fuzz != fm_from.end()) {
                    curr_rg.AddFuzzTo(*subtr_fuzz->second, it.GetStrand());
                }
                const TFuzzMap& fm_to = IsReverse(it.GetStrand()) ? fuzz_to_minus : fuzz_to_plus;
                subtr_fuzz = fm_to.find(curr_rg.GetFrom());
                if (subtr_fuzz != fm_to.end()) {
                    curr_rg.AddFuzzFrom(*subtr_fuzz->second, it.GetStrand());
                }

                if ( have_range  &&  last_id == idh ) {
                    if (x_MergeRanges(last_rg,
                                    last_strand,
                                    curr_rg,
                                    it.GetStrand(),
                                    flags)) {
                        have_range = true;
                        continue;
                    }
                }
                // No merging - push current range, reset last values
                if ( have_range ) {
                    x_PushRange(dst, last_id, last_rg, last_strand);
                }
                last_id = idh;
                last_rg = curr_rg;
                last_strand = it.GetStrand();
                have_range = true;
            }
        }
        else {
            if ( have_range ) {
                bool merged = false;
                if (last_id == idh) {
                    merged = x_MergeRanges(last_rg,
                        last_strand,
                        it_range,
                        it.GetStrand(),
                        flags);
                }
                if ( !merged ) {
                    x_PushRange(dst, last_id, last_rg, last_strand);
                }
            }
            last_id = idh;
            last_rg = it_range;
            last_strand = it.GetStrand();
            have_range = true;
        }
    }
    if ( have_range ) {
        x_PushRange(dst, last_id, last_rg, last_strand);
    }
    if (dst.Which() == CSeq_loc::e_not_set) {
        dst.SetNull();
    }
}


static
void x_SubAndSort(CSeq_loc& dst,
                  const CSeq_loc& minuend,
                  const TFuzzMap& fuzz_from_plus,
                  const TFuzzMap& fuzz_from_minus,
                  const TFuzzMap& fuzz_to_plus,
                  const TFuzzMap& fuzz_to_minus,
                  TIdToRangeColl& rg_coll_plus,
                  TIdToRangeColl& rg_coll_minus,
                  ISynonymMapper& syn_mapper,
                  ILengthGetter& len_getter,
                  CSeq_loc::TOpFlags flags)
{
    bool use_strand = (flags & CSeq_loc::fStrand_Ignore) == 0;

    // Id -> range map for both strands
    unique_ptr<TIdToRangeMap> p_id_map_minus(use_strand ?
        new TIdToRangeMap : 0);
    TIdToRangeMap id_map_plus;
    TIdToRangeMap& id_map_minus = use_strand ?
        *p_id_map_minus.get() : id_map_plus;

    // Prepare default strands
    ENa_strand default_plus = use_strand ?
        eNa_strand_plus : eNa_strand_unknown;
    ENa_strand default_minus = use_strand ?
        eNa_strand_minus : eNa_strand_unknown;

    for (CSeq_loc_CI it(minuend, CSeq_loc_CI::eEmpty_Allow); it; ++it) {
        CSeq_id_Handle_Wrapper idh(syn_mapper.GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        TRangeWithFuzz it_range = TRangeWithFuzz(it);
        if ( it_range.IsWhole() ) {
            it_range.SetOpen(0, len_getter.GetLength(it.GetSeq_id()));
            it_range.ResetFuzzFrom();
            it_range.ResetFuzzTo();
        }
        TRangeColl it_rg_coll(it_range);
        TRanges& rg_map = IsReverse(it.GetStrand()) ?
            id_map_minus[idh] : id_map_plus[idh];
        TIdToRangeColl& rg_coll = IsReverse(it.GetStrand()) ?
            rg_coll_minus : rg_coll_plus;
        TIdToRangeColl::const_iterator id_it = rg_coll.find(idh);
        bool modified = false;
        if (id_it != rg_coll.end()) {
            // Check if there's anything to subtract
            ITERATE(TRangeColl, check_it, id_it->second) {
                if ( it_rg_coll.IntersectingWith(*check_it) ) {
                    it_rg_coll -= id_it->second;
                    modified = true;
                    break;
                }
            }
        }
        if ( modified ) {
            ITERATE(TRangeColl, rg_it, it_rg_coll) {
                TRangeWithFuzz curr_rg(*rg_it);
                // If range start/stop comes from the minuend, preserve strand.
                if (curr_rg.GetFrom() == it_range.GetFrom()) {
                    curr_rg.AddFuzzFrom(it_range);
                }
                if (curr_rg.GetTo() == it_range.GetTo()) {
                    curr_rg.AddFuzzTo(it_range);
                }

                // If range start/stop comes from subtrahend, copy fuzz.
                const TFuzzMap& fm_from = IsReverse(it.GetStrand()) ? fuzz_from_minus : fuzz_from_plus;
                TFuzzMap::const_iterator subtr_fuzz = fm_from.find(curr_rg.GetToOpen());
                if (subtr_fuzz != fm_from.end()) {
                    curr_rg.AddFuzzTo(*subtr_fuzz->second, it.GetStrand());
                }
                const TFuzzMap& fm_to = IsReverse(it.GetStrand()) ? fuzz_to_minus : fuzz_to_plus;
                subtr_fuzz = fm_to.find(curr_rg.GetFrom());
                if (subtr_fuzz != fm_to.end()) {
                    curr_rg.AddFuzzFrom(*subtr_fuzz->second, it.GetStrand());
                }

                rg_map.push_back(curr_rg);
            }
        }
        else {
            rg_map.push_back(it_range);
        }
    }

    x_RangesToSeq_loc(dst, id_map_plus, default_plus, flags);
    if ( use_strand ) {
        x_RangesToSeq_loc(dst, id_map_minus, default_minus, flags);
    }
    if (dst.Which() == CSeq_loc::e_not_set) {
        dst.SetNull();
    }
}


class CDummySynonymMapper : public ISynonymMapper
{
public:
    CDummySynonymMapper(void) {}
    virtual ~CDummySynonymMapper(void) {}

    virtual CSeq_id_Handle GetBestSynonym(const CSeq_id& id)
        {
            return CSeq_id_Handle::GetHandle(id);
        }
};


class CDummyLengthGetter : public ILengthGetter
{
public:
    CDummyLengthGetter(void) {}
    virtual ~CDummyLengthGetter(void) {}

    virtual TSeqPos GetLength(const CSeq_id&)
        {
            return CSeq_loc::TRange::GetWholeToOpen();
        }
};


CRef<CSeq_loc> CSeq_loc::Merge(TOpFlags        flags,
                               ISynonymMapper* syn_mapper) const
{
    unique_ptr<CDummySynonymMapper> p_mapper;
    if ( !syn_mapper ) {
        p_mapper.reset(new CDummySynonymMapper);
        syn_mapper = p_mapper.get();
    }

    CRef<CSeq_loc> ret(new CSeq_loc);
    if ( (flags & CSeq_loc::fMerge_SingleRange) != 0 ) {
        x_SingleRange(*ret, *this, *syn_mapper);
    }
    else if ( (flags & CSeq_loc::fSort) == 0 ) {
        x_MergeNoSort(*ret, *this, flags, *syn_mapper);
    }
    else {
        x_MergeAndSort(*ret, *this, flags, *syn_mapper);
    }
    return ret;
}


CRef<CSeq_loc> CSeq_loc::Add(const CSeq_loc& other,
                             TOpFlags        flags,
                             ISynonymMapper* syn_mapper) const
{
    unique_ptr<CDummySynonymMapper> p_mapper;
    if ( !syn_mapper ) {
        p_mapper.reset(new CDummySynonymMapper);
        syn_mapper = p_mapper.get();
    }

    CRef<CSeq_loc> ret(new CSeq_loc);
    CSeq_loc tmp;
    tmp.SetMix().AddSeqLoc(const_cast<CSeq_loc&>(*this));
    tmp.SetMix().AddSeqLoc(const_cast<CSeq_loc&>(other));
    if ( (flags & CSeq_loc::fMerge_SingleRange) != 0 ) {
        x_SingleRange(*ret, tmp, *syn_mapper);
    }
    else if ( (flags & CSeq_loc::fSort) == 0 ) {
        x_MergeNoSort(*ret, tmp, flags, *syn_mapper);
    }
    else {
        x_MergeAndSort(*ret, tmp, flags, *syn_mapper);
    }
    return ret;
}


CRef<CSeq_loc> CSeq_loc::Subtract(const CSeq_loc& other,
                                  TOpFlags        flags,
                                  ISynonymMapper* syn_mapper,
                                  ILengthGetter*  len_getter) const
{
    unique_ptr<CDummySynonymMapper> p_mapper;
    if ( !syn_mapper ) {
        p_mapper.reset(new CDummySynonymMapper);
        syn_mapper = p_mapper.get();
    }
    unique_ptr<CDummyLengthGetter> p_getter;
    if ( !len_getter ) {
        p_getter.reset(new CDummyLengthGetter);
        len_getter = p_getter.get();
    }

    CRef<CSeq_loc> ret(new CSeq_loc);

    bool use_strand = (flags & CSeq_loc::fStrand_Ignore) == 0;

    // Range collection for each strand
    unique_ptr<TIdToRangeColl> p_rg_coll_minus(use_strand ?
        new TIdToRangeColl : 0);
    TIdToRangeColl rg_coll_plus;
    TIdToRangeColl& rg_coll_minus = use_strand ?
        *p_rg_coll_minus.get() : rg_coll_plus;

    // Collect fuzzes so that they can be copied to the new ranges after subtracting.
    // Note: if there are multiple ranges with the same boundaries only the last
    // fuzz found will be used.
    TFuzzMap fuzz_from_plus, fuzz_from_minus, fuzz_to_plus, fuzz_to_minus;

    // Create range collection(s) for loc2
    for (CSeq_loc_CI it(other); it; ++it) {
        if ( it.IsEmpty() ) {
            continue;
        }
        CSeq_id_Handle_Wrapper idh(syn_mapper->GetBestSynonym(it.GetSeq_id()), it.GetSeq_id());
        TRangeColl& rmap = IsReverse(it.GetStrand()) ?
            rg_coll_minus[idh] : rg_coll_plus[idh];
        rmap += TRangeWithFuzz(it);

        const CInt_fuzz* fuzz = it.GetFuzzFrom();
        if ( fuzz ) {
            if (it.IsSetStrand()  &&  IsReverse(it.GetStrand())) {
                fuzz_from_minus[it.GetRange().GetFrom()].Reset(fuzz);
            }
            else {
                fuzz_from_plus[it.GetRange().GetFrom()].Reset(fuzz);
            }
        }
        fuzz = it.GetFuzzTo();
        if ( fuzz ) {
            if (it.IsSetStrand()  &&  IsReverse(it.GetStrand())) {
                fuzz_to_minus[it.GetRange().GetToOpen()].Reset(fuzz);
            }
            else {
                fuzz_to_plus[it.GetRange().GetToOpen()].Reset(fuzz);
            }
        }
    }

    if ( (flags & CSeq_loc::fMerge_SingleRange) != 0 ) {
        x_SingleRange(*ret,
                      *this,
                      fuzz_from_plus, fuzz_from_minus, fuzz_to_plus, fuzz_to_minus,
                      rg_coll_plus,
                      rg_coll_minus,
                      *syn_mapper,
                      *len_getter);
    }
    else if ( (flags & CSeq_loc::fSort) == 0 ) {
        x_SubNoSort(*ret,
                    *this,
                    fuzz_from_plus, fuzz_from_minus, fuzz_to_plus, fuzz_to_minus,
                    rg_coll_plus,
                    rg_coll_minus,
                    *syn_mapper,
                    *len_getter,
                    flags);
    }
    else {
        x_SubAndSort(*ret,
                     *this,
                     fuzz_from_plus, fuzz_from_minus, fuzz_to_plus, fuzz_to_minus,
                     rg_coll_plus,
                     rg_coll_minus,
                     *syn_mapper,
                     *len_getter,
                     flags);
    }

    return ret;
}


CRef<CSeq_loc> CSeq_loc::Intersect(const CSeq_loc& other,
                                   TOpFlags        flags,
                                   ISynonymMapper* syn_mapper) const
{
    unique_ptr<CDummyLengthGetter> len_getter(new CDummyLengthGetter);
    CRef<CSeq_loc> tmp = Subtract(other,
        // This flag should be used only in the second subtraction
        flags & ~fMerge_SingleRange,
        syn_mapper, len_getter.get());
    return Subtract(*tmp, flags, syn_mapper, len_getter.get());
}


void CSeq_loc::SetStrand(ENa_strand strand)
{
    switch ( Which() ) {
    case e_Int:
        SetInt().SetStrand(strand);
        break;
    case e_Pnt:
        SetPnt().SetStrand(strand);
        break;
    case e_Packed_int:
        SetPacked_int().SetStrand(strand);
        break;
    case e_Packed_pnt:
        SetPacked_pnt().SetStrand(strand);
        break;
    case e_Mix:
        SetMix().SetStrand(strand);
        break;

    default:
        break;
    }
}


void CSeq_loc::ResetStrand(void)
{
    switch ( Which() ) {
    case e_Int:
        SetInt().ResetStrand();
        break;
    case e_Pnt:
        SetPnt().ResetStrand();
        break;
    case e_Packed_int:
        SetPacked_int().ResetStrand();
        break;
    case e_Packed_pnt:
        SetPacked_pnt().ResetStrand();
        break;
    case e_Mix:
        SetMix().ResetStrand();
        break;

    default:
        break;
    }
}


END_objects_SCOPE // namespace ncbi::objects::
END_NCBI_SCOPE

#undef NCBI_USE_ERRCODE_X
